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INTRODUCTION 

As  a  body  of  revolution  is  exposed  to  an  increasing  angle  of 
attack,  the  flow  on  its  leeward  side  separates  and  rolls  up  to  form 
vortices.  Under  these  conditions.,  normal  force  and  pitching  moment 
can  no  longer  be  predicted  by  slender  body  theory.  For  relatively 
small  angles  of  attack  (a  <  tan“l (4 .5/£) )  two  symmetric  vortices 
form  as  shown  in  Figure  la.  At  larger  angles  of  attack  wake  vortices 
grow  unsymmetrically  and  break  away  or  are  shed  from  the  leeward 
side,  as  illustrated  in  Figure  lb.  The  resulting  flow’-field 
asymmetries  induce  a  yawing  force  and  moment. 

The  analysis  of  bodies  of  revolution  at  angles  of  attack  has 
been  approached  by  many  workers  with  varying  degrees  of  rigor. 

Simple  engineering  calculations  can  be  used  to  determine  normal 
forces  and  pitching  moments  (Refs.  (l)-(3)),  but  extensive  numerical 
computations  must  be  employed  if  detailed  flow-field  information  is 
desired  (Refs.  (4),  (5)).  The  objective  of  this  report  is  to  predict 

the  yawing  force  component.  Since  this  force  arises  from  asymmetries 
in  the  vortex  wake,  the  appropriate  level  of  analysis  must  be 
sufficiently  rigorous  to  provide  an  approximate  wake  description. 
Engineering  level  calculations  do  not  generally  provide  this 
description  and,  therefore,  are  not  amenable  to  these  ends.  While 
detail  flow-field  calculations  can,  in  principle,  achieve  this 
purpose,  they  have  only  been  applied  up  to  very  moderate  angles  of 
attack  (10-15  degrees) ,  too  low  for  yawing  forces  to  occur.  Semi- 
empirical  methods,  such  as  those  by  Bryson  (Ref.  (6)) ,  Schindel 
(Ref.  (7)),  and  Angelucci  (Ref.  (8)),  provide  an  approximate  wake 
description,  but  assume  a  symmetric  vortex  wake  and  are  not  capable 
of  predicting  yawing  forces.  Also,  normal  forces  and  pitching 
moments  can  only  be  calculated  up  to  moderate  angles  of  attack 
(20-30  degrees) . 

The  theory  developed  in  this  paper  extends  the  mode],  nri  rri-np l  ju, 
proposed  by  Bryson  (Ref.  (6))  by  removing  the  symmetric  wake 
assumption.  This  allows  the  side  forces  to  be  predicted,  and 
increases  angle  of  attack  for  which  the  theory  is  valid  to  about 
50  degrees.  In  the  current  theory,  the  flow  field  is  approximated 
using  the  impulsive  flow  analogy.  This  analogy  assumes  that  the 
crossflow  field  is  swept  down  the  length  of  the  body  at  the  rate 
Ucosa.  At  each  axial  station  the  flow  field  is  taken  to  be 
analogous  to  flow  about  a  cylinder  in  crossflow,  whose  radius  is 
equal  to  the  body  radius  at  that  axial  location.  The  developing 
crossflow  field  is  thus  similar  to  flow  about  an  impulsively 
started  cylinder  in  crossflow.  The  development  of  the  viscous  flow 
field  is  modeled  by  superimposing  point  vortices  on  the  potential 
flow  solution  for  flow  about  a  cylinder.  Analysis  is  limited  to 
ogive  cylinders  with  sharp  or  blunt  nose  tips. 

Using  experimentaldata  for  a  body  with  a  laminar  boundary 
layer  (2.5  <  Re/ft  x  10  <  5.6)  and  covering  a  range  in  Mach  number 

from  0.5  to  1.1,  free  parameters  inherent  in  the  model  are  determined. 
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Calculated  force  and  moment  values  are  compared  to  experiment  and 
computed  vortex  paths  are  compared  to  those  determined  from  schlieren 
photographs.  This  type  of  approach  allots  conclusions  to  be  drawn 
concerning  the  validity  of  the  current  problem  formulation. 
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local  drag  coefficient,  j 

sin  a 


M 


pitching  moment  coefficient,  - ■L— 

model  nose  Dq  nr 


2,  referenced  to 


M 


yawing  moment  coefficient, 


normal  force  coefficient, 


Dq^r 
F» 


N 


CT  TTr 

~oo 


JL 

yawing  force  coefficient,  — 

q  irr 4 

T.00 

diameter  of  body  of  revolution 
normal  force 


yawing  force 

number  of  vortices  in  wake 
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body  fineness  ratio 
nose  fineness  ratio 

Mach  number 

crossflow  Mach  number,  Msina 
pitching  moment 
yawing  moment 

frequency  for  shedding  vortices  of  like  sign 


2 


NOLTR  73-209 


free-stream  dynamic  pressure 


R  nose  tip  radius  of  curvature 
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r^  maximum  body  radius 
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difference  in  the  radial  location  of  the  two  initial 
vortices  at  which  the  first  vortex  will  be  shed. 

the  distance  from  the  separation  point,  measured  along  the 
feed  sheet  at  which  nascent  vortices  are  introduced 

Reynolds  number 

perturbation  in  the  radial  location  of  the  first  two 
vortices 
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Strouhal  number 
time 

free-stream  velocity 
crossflow  plane  velocity  at  t,, 


d<HC) 


dimensionless  crossflow  plane  velocity,  W/Usina 

dimensionless  crossflow  plane  velocity  at  the  location 
of  the  jth  vortex.  The  influence  of  the  jth  vortex 
is  neglected. 

distance  along  missile  body  (see  Fig.  1) 

x  distance  between  points  where  successive  vortices  are 
shed. 

crossflow  plane  coordinates  (see  Fig.  2) 
y  and  z  coordinates  of  the  jth  vortex 
angle  of  attack 
vortex  circulation 

.  dx 

rate  of  radius  change  as  a  function  of  x,  ^ 
complex  crossflow  plane  coordinate,  y  +  iz 
crossflow  plane  separation  points  (see  Fig.  2) 
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crossflow  plane  separation  angles  (see  Fig.  2) 

•  ,  P 

dimensionless  circulation,  — i _ 

2irrsinaU 

dimensionless  circulation  of  newly  introduced  vortice 
initial  vortex  strength 
crossflow  velocity  potential 
nose  roll  angle 

dimensionless  crossflow  velocity  potential,  4>/Usina 


LITERATURE  SURVEY 

The  complete  analysis  of  a  body  of  revolution  at  angle  of  attack 
is  a  three-dimensional  problem  involving  the  solutions  to  the 
Navier-Stokes  equations.  This  problem  may  be  approached  numerically 
in  its  entirety  requiring  an  extensive  computational  effort.  For 
relatively  small  angles  of  attack  (12  degrees) ,  numerical  solutions 
have  been  obtained  for  cones  using  an  approximate  set  of  Navier- 
Stokes  equations  (Ref.  (4)). 

The  device  most  commonly  used  to  simplify  the  problem  formulation 
is  the  impulsive  flow  analogy.  This  reduces  the  complete  problem  to 
a  two-dimensional  time  varying  analysis  of  a  viscous  fluid.  By  this 
analogy,  the  crossflow  field  is  equivalent  to  that  about  an  impulsively 
started,  expanding  cylinder,  where  t  =  x/Ucos ,  The  expanding 
cylinder  represents  the  increase  in  cross-sectional  radius  which 
occurs  along  the  length  of  the  body.  Often  the  additional  assumption 
is  made  that  the  flow  about  an  impulsively  started  cylinder  is 
incompressible.  Throughout  the  remainder  of  this  report  the  term 
impulsive  flow  analogy  will  be  assumed  to  imply  this  assumption. 

On  the  most  rigorous  level,  the  impulsive  flow  analogy  allows 
the  crossflow  drag  to  be  deduced  from  a  numerical  solution  to  the 
Navier-Stokes  equations  for  flow  over  an  impulsively  started  cylinder 
in  crossflow.  Several  examples  of  this  type  of  procedure  are 
available  for  cylinders  of  constant  radius  (Refs.  (9),  (10)).  For 

direct  application  to  a  body  at  angle  of  attack  Walitt  and  Trulio 
(Ref.  (5))  have  carried  out  this  analysis  with  a  cylinder  of  varying 
radius  in  compressible  flow.  However,  the  computational  requirements 
are  prohibitive  and  solutions  have  only  been  obtained  for  angles  of 
attack  of  15  degrees. 

A  less  rigorous  analysis  of  impulsively  started  flow  over  a 
cylinder  involves  superposition  of  a  large  number  of  point  vortices 
on  the  potential  solution  for  incompressible  flow  about  a  cylinder. 
Typically,  several  hundred  discrete  vortices  may  be  used.  After 
formation,  each  point  vortex  is  assumed  to  drift  with  the  local 
fluid  velocity.  The  schemes  for  predicting  rate  of  vorticity 
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formation  differ  extensively.  Although  many  such  studies  have  been 
carried  out,  most  are  restricted  to  bodies  of  fixed  cross-sectional 
dimensions  (Refs,  (11),  (12)  and  (13)),  However,  Angelucci  (Ref.  (8)) 
has  made  a  direct  application  to  bodies  at  angle  of  attack,  assuming 
a  symmetric  vortex  wake.  Good  agreement  is  obtained  for  normal 
forces  and  moments. 

A  somewhat  simpler  method  uses  a  single  vortex  to  describe  each 
wake  vortex.  In  a  procedure  originally  proposed  by  Bryson  (Ref.  (6)) , 
vortex  motion  is  determined  by  a  balance  of  forces  on  a  vortex  and 
its  feeding  sheet.  The  model  was  originally  limited  to  circular 
bodies  with  a  wake  pattern  containing  two  symmetric  vortices.  This 
approach  has  been  extended  by  Schindel  (Ref.  (7))  to  bodies  of 
elliptic  cross  section  using  conformal  mapping  techniques.  In 
addition,  changes  in  empirical  parameters  allow  laminar  and  turbulent 
boundary  effects  to  be  accounted  for.  Results  compare  well  with 
experimental  data  for  normal  forces.  However,  agreement  of  the 
pitching -moment  calculation  with  experimental  data  is  not  as  good. 

Range  of  validity  for  the  results  is  a  <  tan  1(5/£).  Davis  (Ref. 

(14))  has  extended  Bryson's  methods  to  a  cylinder  in  crossflow  with 
an  asymmetric  wake.  But  this  formulation  is  valid  only  for  cylinders 
of  fixed  radius  and,  therefore,  cannot  be  directly  applied  to  a  body 
at  angle  of  attack. 

The  most  approximate  theories  available  consider  forces  on  a 
body  at  angle  of  attack  to  be  made  up  of  a  viscous  and  an  inviscid 
contribution  which  are  additive.  Initial  attempts  along  these  lines 
were  by  H.  J.  Allen  (Ref.  (1)).  The  inviscid  force  expression 
applied  is  similar  to  that  obtained  from  slender  body  theory.  The 
viscous  contribution  at  each  axial  station  is  computed  by  applying 
the  steady-state  drag  force  for  a  cylinder  in  crossflow.  This  report 
was  followed  by  that  of  Kelly  (Ref.  (2  ))  who  attempted  to  improve 
the  method  by  using  drag  data  of  an  impulsively  started  cylinder  in 
crossflow.  His  results  were  somewhat  better  than  those  of  Allen,  giving 
a  better  prediction  of  the  pitching  moment.  Kelly's  method  has  been 
refined  by  others,  most  recently  by  Thomson  (Ref.  (15))  to  include 
angles  of  attack  up  to  90  degrees.  In  this  work  advantage  is  taken 
of  impulsive  flow  data  obtained  by  Sarpkaya  (Ref.  (16)) .  Correction 
tables  allow  Mach  number,  end  effects,  nose  shapes,  and  the  transition 
from  laminar  to  turbulent  flow  to  be  taken  into  consideration.  In 
addition,  an  experimental  wake  description  is  developed  which  allows 
side  forces  to  be  estimated.  In  a  slightly  different  vein,  L.  H. 
Jorgensen  (Ref.  (3))  in  an  extension  of  Allen's  work,  presents  a 
method  of  estimating  normal  forces  on  bodies  at  an  angle  of  attack 
varying  from  0  to  180  degrees  and  Mach  numbers  varying  from  subsonic 
to  hypersonic. 
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BRYSON'S  METHOD  FOR  IMPULSIVELY  STARTED  CYLINDERS 
IN  INCOMPRESSIBLE  FLOW 

As  pointed  out  above,  our  analysis  is  essentially  that  of  an 
impulsively  started  cylinder  in  incompressible,  viscous  flow.  Such 
flow  can  be  simulated  by  superimposing  point  vortices  on  the  potential 
solution  for  flow  about  a  cylinder.  In  the  present  approach,  a 
single  point  vortex  is  used  to  simulate  each  wake  vortex.  To  satisfy 
the  boundary  conditions  on  a  cylinder,  each  vortex  in  the  wake  is 
required  to  have  a  corresponding  image  vortex.  The  resulting  con¬ 
figuration  is  shown  in  Figure  2  and  the  complex  velocity  potential 
is  given  as  follows : 
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The  term , 


is  the  velocity  potential  solution  for 


flow  over  a  cylinder.  Each  of  the  logarithic  terms  within  the  sum¬ 
mation  represent  the  velocity  potential  for  a  vortex  and  its  image 
within  the  cylinder.  The  remaining  term,  which  is  the  expression 
for  a  source,  takes  into  account  the  changing  cylinder  radius. 


The  equations  governing  the  motion  of  the  vortices  remain  to  be 
determined.  Each  logarithmic  term  within  the  summation  in  Equation 
(1)  is  a  double  valued  function  and  has  a  branch  line  which  connects 
the  described  vortex  to  its  image.  The  portion  of  the  branch  line  of 
physical  importance  is  the  part  which  crosses  the  flow  field.  It  is 
considered  to  extend  from  the  vortex  location  to  the  nearer  of  the 
points,  and  T,  ,  which  are  on  the  cylinder  surface  (see  Fig.  2). 

1  U2 


These  two  points  correspond  to  the  position  where  flow  is  assumed  to 
separate  from  the  cylinder  surface.  The  branch  line  of  a  growing 
vortex  (dT/dt  Y  0)  is  known  as  a  feeding  sheet.  Across  it,  a 

pressure  jump  of  occurs,  which  can  be  determined  from  the 

momentum  conservation  equation.  The  total  force  on  the  feeding  sheet 

is  ip  “  Cq  force  on  the  feeding  sheet  is  balanced  by  a 

k  a?k 

lift  force  on  the  vortex  of  -ipr  -  U^.)  given  by  the  Kutta- 


Joukowski  theorem.  Here  (■ 


dt 


-  U,  )  is  the  vortex  velocity  with 
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respect  to  the  local  fluid  velocity.  Thus,  the  motion  of  each  vortex 
is  governed  by  the  equation 


££■  +  (f,  -  s.j ) 

cU  n.  eft 


*? 


= 


(2) 


The  vortex  circulations  are  determined  by  applying  the  Kutta 
conditions  at  the  crossflow  separation  points: 


W(~Sa, )  =  W(  % oO  =  o 


(3) 


This  allows  the  circulation  of  two  different  vortices  at  a  time  to 
be  specified.  This  does  not  preclude  the  existence  of  other  vortices 
in  the  flow  whose  circulations  are  known. 

Left  undetermined  by  the  Bryson  method  are  the  crossflow 
separation  points,  and  Coo'  and  boundary  conditions  in  t*  The 

boundary  conditions  t  include  the  points  where 

vortices  are  introduced  and  shed.  Also,  the  location  of  vortices 
in  the  crossflow  plane  and  their  strengths  must  be  specified  when 
they  are  introduced . 

In  applying  the  Bryson  method  the  boundary  conditions  in  t 
must  be  specified  in  such  a  manner  that  the  resulting  wake  conforms 
to  observations.  When  the  wake  behind  the  cylinder  is  assumed  to 
become  asymmetric,  the  following  type  of  procedure  appears  to  be 
required.  Initially,  two  vortices  are  introduced.  Using  the 
governing  equations  the  growth  and  the  motion  of  these  vortices 
are  traced  numerically.  At  some  point,  to  be  empirically  deter¬ 
mined,  the  circulation  of  one  of  these  vortices  is  fixed  and  this 
vortex  is  considered  to  be  shed.  A  new  vortex  is  then  introduced 
on  the  same  side  of  the  cylinder  that  the  shed  vortex  originally 
emanated  from.  This  shedding  process  is  repeated  at  later  points 
in  time,  vortices  being  shed  from  alternate  sides  of  the  cylinder. 


APPLICATION  OF  THE  BRYSON  METHOD  TO  A  MISSILE  AT  ANGLE  OF  ATTACK 

The  impulsive  flow  analogy  implies  the  relations  x=Ucosat  and 
W(ro)  =  Usina.  Using  these  expressions  and  the  definitions 
f  =  2irUr  sinaA  and  Wj  =  WjUsina,  Equation  (2)  becomes 


-f-  (  Tj  -  )  ciA  _  tcMrtu;  -  ( )  dr 

A  _j'  dP  3  r  dp 

—  tam  c<  <jl >  .j  j  h  *  •  ■  K 

d-* 


shed 


growing 

vortices 


vortices 
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The  Kutta  condition  for  a  cylinder  whose  radius  may  be 
expanding  becomes : 


1 

1«,)  _ 

UJ  (  lot  )  - 

r 

dr  J_ 

d  3  To  | 

1 

d  4  « 

(  lo 1  )  . 

r 

dr  ,j_t 

LA 

'+O^VveC 

d2:  f 

Substituting  Equation  (1)  into  Equations  (5) ,  it  becomes  evident 
that  Xq  an<^  ^2  are  determined  by  solving  two  simultaneous  equations. 
It  is  also  evident  in  this  formulation  that  satisfaction  of  the 
inviscid  boundary  condition  on  the  cylinder  is  sufficient  to 
ensure  that  both  Xq  and  X2  are  real. 

In  solving  Equations  (4)  subject  to  conditions  (5),  Xq  and  X2 
have  been  eliminated  from  (4)  reducing  the  number  of  equations 
involved.  The  actual  derivation  is  shown  in  Appendix  A. 

The  formulated  method  is  applied  in  the  following  manner. 

Near  the  nose  of  the  missile  two  nascent  vortices  are  introduced. 

The  x  coordinate  of  introduction  is  determined  by  satisfying  the 
relation  —  sin0  tana  =  tan6 .  This  relationship  moves  the  point 
of  vortex^introduction  toward  the  nose  of  the  model  with  increasing 
angle  of  attack.  In  cases  where  this  relation  places  the  starting 
point  in  front  of  the  nose,  the  problem  is  automatically  started 
at  the  nose.  This  relationship  was  derived  by  Bryson  (Ref.  (6)) 
by  considering  a  cone.  The  slope  6  represents  the  half  angle  of 
the  largest  cone  for  which  vortex  separation  occurs.  The  two 
vortices  are  each  assumed  to  be  located  along  the  directions 
60  +  1T/3  and  0q2  -  tt/3  from  their  respective  separation  points 

6q^  and  <5q2-  This  criterion  is  also  due  to  Bryson  and  was  deter¬ 
mined  by  the  analysis  of  the  trajectory  of  a  vortex  near  the 
crossflow  separation  point. 

Asymmetry  is  introduced  at  this  point  by  slightly  perturbing 
the  radial  location* of  the  two  vortices.  The  initial  strengths  of 
the  two  vortices,  as  well  as  the  magnitude  of  the  perturbation 
applied,  are  free  parameters  which  are  empirically  determined. 

They  will  be  discussed  in  more  detail  later. 


As  the  numerical  solution  to  the  governing  equations  progresses, 
the  asymmetry  in  the  problem  grows.  When  the  difference  in  radial 
location  of  the  two  vortices  reaches  a  specified  value,  the 
outermost  vortex  is  shed  and  a  new  vortex  is  introduced  on  the  same 

*the  radial  location  of  the  jth  vortex  is  defined  to  be 
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side  of  the  body.  In  order  to  introduce  a  new  vortex  whose  circula¬ 
tion  is  other  than  zero,  the  existing  vortex  pattern  must  be 
perturbed.  This  is  accomplished  by  reducing  the  circulation  of 
the  shed  vortex  a  few  percent.  It  can  be  shown  that  this  procedure 
will  always  lead  to  a  new  vortex  with  a  circulation  of  the  correct 
sign.  The  new  vortex  is  assumed  to  be  located  along  the  feeding 
sheet  of  the  vortex  to  be  shed.  Its  initial  strength  is  uniquely 
determined  by  its  crossflow  field  location  and  the  decrease  in 
circulation  ascribed  to  the  shed  vortex. 


Subsequent  vortices  are  shed  at  intervals  specified  by  the 
Strouhal  frequency.  Initial  conditions  for  succeeding  vortices 
are  identical  to  those  just  described. 


Application  of  the  method  described  above  requires  specification 
of  the  following  free  constants : 

1.  Initial  vortex  strength  (Aj) 

2.  Initial  radial  perturbation  (Ar  ) 

P 

3.  The  difference  in  radial  location  of  the  two  initial  vortices 
at  which  the  first  vortex  is  shed  (Ar) 


4.  The  distance  from  the  separation  point,  measured  along  the 
feeding  sheet  of  the  shed  vortex,  at  which  a  nascent  vortex  is 
introduced  (rn) 

5.  Circulation  reduction  of  the  shed  vortex  (AAn) 

6.  Crossflow  separation  angles  (0q^  and  Gg^) 

7.  Strouhal  number  (S) 


The  Strouhal  number  is  determined  experimentally  using  schlieren 
photographs  of  a  tangent  ogive  cylinder  at  angle  of  attack.  Although 
measured  values  showed  a  large  degree  of  fluctuation,  an  average 
value  of  .17  was  established. 

The  crossflow  separation  angle  can  also  be  determined  in  the 
wind  tunnel,  however,  use  of  measured  values  leads  to  unstable  numerical 
results.  This  may  be  due  to  the  approximate  nature  of  the  flow 
field  model.  Simulation  of  a  wake  vortex  using  a  single  point 
vortex  may  not  be  giving  an  accurate  flow  field  description  in  the 
vicinity  of  the  separation  points.  To  circumvent  this  problem,  the 
lead  of  Bryson  is  followed  and  the  separation  angle  is  set  near 
45  degrees. 

There  is  no  direct  way  of  readily  determining  the  remaining 
parameters.  However,  it  should  be  pointed  out  that  the  numerical 
solution  is  much  less  sensitive  to  these  remaining  parameters. 

Their  actual  values  were  established  by  a  comparison  with  experi¬ 
mental  measurements  of  normal  forces  and  schlieren  pictures  of  the 
leeward  flow  field. 
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NORMAL  AND  YAWING  FORCES 

The  forces  acting  on  any  section  of  a  body  can  be  determined 
by  enclosing  it  in  a  cylindrical  control  volume  and  balancing  the 
momentum  flux  through  the  volume.  Such  a  procedure  is  adequately 
discussed  elsewhere  (Ref.  (17))  and  results  in  the  following 
expression: 

k 

Cfj  f  I  Cy  -  Z  ^  £  J*  £  /  f  2  -  JZ  J  J 

The  first  term  in  the  above  equation  reduces  to  the  slender 
body  theory  result  of  2a  at  small  angles  of  attack.  The  second 
term  represents  the  crossflow  drag  due  to  the  presence  of  the 
vortices,  or  the  viscous  contribution. 


EVALUATION  OF  THE  PROBLEM  FORMULATION 

The  validity  of  the  current  formulation  hinges  on  several 
points  which  are  difficult  to  resolve  in  a  conclusive  manner. 

First  and  foremost,  is  the  validity  of  the  impulsive  flow  analogy. 
This  analogy  cannot  be  justified  in  its  entirety.  However,  for 
inviscid  flow  about  a  slender  body  at  small  angle  of  attack,  it  is 
a  consequence  of  slender  body  theory.  At  low  Mach  numbers  where 
the  flow  is  approximately  incompressible,  the  restriction  on  small 
angle  of  attack  can  be  relaxed.  This  does  not,  of  course,  shed  any 
light  on  the  validity  of  the  impulse  flow  analogy  for  fluid  flows 
in  which  viscous  effects  are  important.  All  that  can  really  be 
said  in  these  cases  is  that  the  observed  development  of  a  vortex 
wake  on  the  leeward  side  of  a  body  at  angle  of  attack  bears  at 
least  a  qualitative  relation  to  the  wake  which  develops  behind  a 
cylinder  in  crossflow  (Ref.  (18)). 

Another  point  of  question  in  the  current  problem  formulation 
is  the  practice  of  simulating  viscous  flow  by  superimposing  point 
vortices  on  inviscid,  potential  flow  solutions.  Even  though  this 
is  an  ad-hoc  problem  formulation,  previous  success  in  studies  for 
cylinders  in  crossflow  using  a  large  number  of  point  vortices 
indicate  that  this  manner  of  approach  provides  an  accurate  wake 
description  -  in  one  case,  good  enough  to  predict  the  shedding 
frequency  (Ref.  (11)),  The  use  of  a  single-point  vortex  to  describe 
each  wake  vortex,  as  is  done  in  this  report,  is  more  questionable. 
Such  an  approach  is  only  capable  of  describing  the  gross  flow-field 
features.  Pressure  calculations  around  the  circumference  of  a 
cylinder  using  this  type  of  description  show  an  unrealistic  peak 
near  each  point  vortex.  Despite  this  shortcoming,  previous  studies 
using  this  approach  indicate  that  the  integrated  circumferential 
pressure  distribution  gives  a  good  approximation  of  the  normal 
forces  at  low  angles  of  attack  (Refs,  (6),  (7)), 
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The  present  approach  is  thus  the  extension  of  other  studies 
which  have  proved  to  be  accurate  at  lower  angles  of  attack.  Free 
parameters  available  in  the  current  model  allow  adjustments  to  be 
made  for  effects  which  are  not  directly  accounted  for. 


NUMERICAL  CALCULATIONS 

The  governing  equations  are  integrated  using  an  Adams-Moulton 
predictor-corrector  scheme  which  adjusts  step  size  based  on 
estimated  error.  This  routine  is  started  using  the  Runge-Kutta 
technique.  The  adjustable  step  size  feature  is  extremely  important 
in  reducing  running  time.  Small  step  sizes  are  required  near  the 
model  nose  where  the  local  radius  is  small  and  at  points  at  which 
new  vortices  are  introduced.  A  listing  of  the  program  is  provided  in 
Appendix  B. 

Correct  use  of  the  outlined  method  requires  the  determination 
of  a  set  of  free  parameter  values  that  not  only  give  a  correct 
qualitative  description  of  the  wake,  but  also  produce  realistic 
forces  and  moments.  There  does  not  appear  to  be  a  unique  set  of 
parameters  that  best  satisfy  these  requirements.  However,  certain 
sets  of  parameter  values  produce  unacceptable  results  in  that  the 
computed  crossflow  wake  does  not  develop  into  a  Karman  vortex  street 
type  of  pattern.  In  some  instances  newly  introduced  vortices  may  move 
with  diminishing  circulation  toward  the  cylinder  surface.  In  any 
event,  unless  an  orderly  wake  is  formed,  point  vortices  eventually 
approach  one  another  resulting  in  numerical  instability. 

Numerical  results  were  initially  calculated  for  an  impulsively 
started  cylinder  in  crossflow.  Suitable  results  were  obtained  using 
the  set  of  parameters  shown  in  Table  1.  As  indicated  in  this  table, 
calculations  were  carried  out  assuming  a  symmetric  wake  containing 
two  vortices,  an  unsymmetric  wake  still  containing  only  a  vortex 
pair,  and  an  asymmetric  wake  in  which  vortices  were  periodically 
shed.  In  Figure  3  the  calculated  drag  curves  are  plotted  as  a 
function  of  dimensionless  time.  Also  shown  in  this  graph  is  the 
experimental  drag  curve  determined  on  an  impulsively  started  cylinder 
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FREE  PARAMETER  VALUES 

_ _  Cylinde r  in  Crossflow  _ _ _ _ _ 

Asymmetric  Vortex  Pair  Asymmetric  wake  with 

Parameter  Symmetric  Vortex  Pair _ no  shedding  _ vortex  shedding 


Missile  at 
angle  of  attack 
Mach  no.  taken 
into  account 
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in  crossflow  by  Sarpkaya  (Ref.  (16)).  Only  the  calculation  which 
provides  for  vortex  shedding  is  comparable  to  experimental  results 
over  wide  ranges  in  dimensionless  time.  The  accompanying  wake 
development  is  illustrated  in  Figure  4. 

In  our  initial  attempt  to  apply  this  method  to  a  body  at  angle 
of  attack,  the  same  set  of  free  parameter  values  used  for  the 
impulsively  started  cylinder  was  assumed.  Results  were  compared 
to  experimental  tests  on  a  tangent  ogive  cylinder  (i '=  12  and  =  4) 
These  tests  are  further  described  in  Appendix  C.  The  only 
adjustment  to  the  set  of  free  parameter  values  concerned  the  cross- 
flow  separation  angles.  As  can  be  seen  in  Table  1  for  the  cylinder 
in  crossflow,  the  separation  angles  are  decreased  at  the  dimension¬ 
less  time  _  =  10.2.  In  application  to  the  ogive  cylinder  at  angle 
of  attack, rthe  point  in  the  computation  where  this  decrease  occurs 
must  be  adjusted,  since  the  change  in  radial  dimensions  at  the 
missile  nose  disturbs  the  time  sequence  of  the  flow  development. 

It  is  found  that  the  point  at  which  =  10.2  occurs  two-thirds 
of  the  distance  in  time  between  the  Joint  where  the  first  and  second 
vortices  are  shed.  This  criterion  is  used  to  determine  where 
0q  and  8q  will  be  decreased.  As  suggested  by  Schindel  (Ref.  (7)), 

the  smooth  transition  from  one  separation  angle  to  another  is 
achieved  by  solving  the  equation: 

d  Qq  j  _  /7  tana.  (  ~  Qoj  )  ^  2  (6) 

d  ji  r 

and  subject  to  the  boundary  condition:  x  ->•  00  as  0q_.  -*  48°.  In 

the  special  case  of  a  cylinder  in  crossflow,  r  is  constant  and 
the  solution  to  16  reduces  to  that  shown  in  Table  1  for  a  cylinder 
in  crossflow. 

The  calculated  forces  and  pitching  moments  are  compared  to 
experimental  results  in  Figures  5  and  6.  In  each  of  these  figures 
experimental  data  for  a  range  in  Mach  number  from  .5  to  1.1  are 
shown.  The  computed  normal  force  coefficient  (Fig.  5)  compares 
well  to  the  experimental  data  being  closest  to  the  curve  generated 
at  Mach  .9.  The  calculated  pitching  moment  coefficient  (F'ig.  6) 
is  slightly  high,  matching  best  with  the  Mach  1.1  curve.  The 
predicted  yawing  force  coefficient,  shown  in  Fig.  7,  is  somewhat 
larger  than  those  measured  and  is  closest  to  the  Mach  .5  curve. 

The  current  theory  is  independent  of  Mach  number.  However, 
as  can  be  seen  from  experimental  data  in  Figures  5,  6,  and  7, 
force  and  moment  data  is  Mach  number  dependent.  In  order  to 
include  a  Mach  number  effect  in  the  theory,  the  set  of  free  con¬ 
stants  is  assumed  to  be  a  weak  function  of  Mach  number.  The 
relation  between  Mach  number  and  free  constant  value  is  determined 
empirically  by  adjusting  the  set  of  constants  to  give  the  measured 
normal  force  curve  at  each  Mach  number.  In  addition,  the  sets  of 
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free  parameter  values  are  modified  somewhat  to  obtain  better  corres¬ 
pondence  between  observed  and  calculated  vortex  trails.  The  final 
set  of  constants  is  shown  in  the  right-hand  column  of  Table  1. 

These  constant  values  yield  stable  vortex  streets  up  to  angles  of 
attack  near  50  degrees.  The  upper  limit  in  each  case  is  dependent 
on  the  actual  nose  configuration.  Generally,  the  upper  limit  on 
angles  of  attack  increases  with  decreasing  nose  fineness. 

Computations  using  the  Mach  number  dependent  set  of  constants 
are  compared  to  experiment  in  Figures  8  through  19.  The  close 
agreement  between  experimental  and  computed  normal  force  coefficients 
are  shown  in  Figure  8  while  the  pitching  moment  results  are  compared 
in  Figure  9.  The  computed  yawing  force  coefficient  is  compared  to 
experimental  measurements  at  Mach  .5,  .7,  .9,  and  1.1,  in  Figures 

10,  11,  12,  and  13,  respectively.  In  each  of  these  graphs  four 
different  yawing  force  coefficient  curves  are  shown.  They  were 
generated  by  rolling  the  model  nose.  From  these  graphs  it  is 
evident  that  the  yawing  force  coefficient  is  reasonably  well  pre¬ 
dicted  for  Mach  numbers  less  than  .9.  For  larger  Mach  numbers  the 
measured  magnitudes  of  the  yaw  force  coefficients  decreases.  This 
trend  is  not  followed  by  the  computed  results .  At  the  highest 
Mach  number  tested,  1.1,  the  yawing  force  coefficient  is  over¬ 
estimated  by  a  factor  of  four. 

In  Figures  14  through  19,  comparisons  between  the  computed 
and  observed  vortex  wakes  are  shown  for  angles  of  attack  varying 
from  20  degrees  to  45  degrees  and  Mach  numbers  of  .5  and  .9.  It 
is  evident  from  these  pictures  that  the  wake  is  qualitatively  well 
reproduced.  More  consistent  agreement  is  obtained  at  lower  Mach 
numbers.  At  an  angle  of  attack  of  45  degrees  and  a  Mach  number 
of  .9,  the  wake  width  is  clearly  over-estimated. 

In  order  to  assess  the  general  applicability  of  this  model, 
a  comparison  is  made  between  calculated  and  experimental  results 
for  differently  shaped  ogive  cylinders.  Results  at  Mach  .5  are 
shown  in  Figures  20  through  22.  Good  agreement  is  obtained  for 
normal  force.  Although  the  pitching  moment  coefficient  comparison 
is  not  as  good  as  the  normal  force  comparison,  the  experimental 
trend  is  followed.  Figure  22  reveals  that  the  decreasing  yawing 
force  coefficient  observed  on  blunter  bodies  is  not  predicted  by 
the  theory . 

Comparison  between  experimental  and  calculated  results  on 
differently  shaped  bodies  at  Mach  1.1  (Figures  23  through  25)  shows 
that  general  trends  found  in  experimental  data  for  normal  force 
and  pitching  moments  are  followed.  However,  agreement  is  not  as 
close  as  in  the  case  of  lower  Mach  numbers.  Also  the  side  forces 
(Fig.  26)  are  consistently  over-estimated  and  do  not  follow  experi¬ 
mental  trends. 
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SUMMARY  AND  CONCLUSIONS 

The  flow  field  surrounding  an  ogive  cylinder  at  high  angles  of 
attack  is  simulated  using  an  extension  of  the  method  initially 
proposed  by  Bryson  (Ref.  (6)).  This  method  employs  the  impulsive 
flow  analogy  and  models  each  wake  vortex  using  a  single  point 
vortex.  Free  parameters  inherent  in  the  problem  formulation  are 
initially  determined  by  comparison  to  drag  force  coefficient  data 
on  an  impulsively  started  cylinder  in  crossflow. 

Final  parameter  values  are  obtained  by  matching  the  normal 
force  measured  on  a  tangent  ogive  cylinder  ( Jin  =  4.0,  i  =  12)  over 
a  range  in  Mach  number  of  .5  to  1.1.  These  parameter  values  are 
made  weakly  Mach  number  dependent  to  reflect  the  experimental  data. 
The  resulting  pitching  moments  compare  well  in  trend  with  experiment. 
The  calculated  yaw  force  coefficient  is  of  the  correct  magnitude 
at  Mach  numbers  less  than  .9,  but  is  too  large  at  higher  Mach 
numbers . 

At  a  Mach  number  of  .5  calculated  and  experimental  results 
are  compared  for  ogive  cylinders  with  different  nose  fineness 
ratios.  Predicted  normal  force  and  pitching  moment  coefficients 
are  in  good  agreement  with  experimental  results.  On  ogive  cylinders 
with  small  slenderness  ratios  (£n  <  3)  yawing  force  calculations 
over-estimated  the  experimental  results. 

The  same  type  of  comparison  is  carried  out  at  a  Mach  number 
of  1.1.  Agreement  between  calculated  normal  force  and  pitching 
moment  coefficients  was  not  as  good  as  at  Mach  .5.  The  yawing 
force  coefficient  was  consistently  lower  than  predicted. 

The  current  results  indicate  that  on  slender  bodies  and  at 
low  Mach  numbers,  the  impulsive  flow  analogy  combined  with  the 
Bryson  method  can  be  used  to  predict  normal  force,  pitching  moment 
and  yawing  force.  For  nose  fineness  ratios  less  than  4  and  at 
Mach  numbers  near  unity,  accuracy  deteriorates  markedly  and  side 
forces  are  over-estimated. 

Apparently  the  separation  parameters  are  sensitive  to  Mach 
number  in  the  transonic  range.  A  more  thorough  investigation  of 
this  effect  would  be  required  to  improve  the  accuracy  of  side 
force  predictions  in  this  regime. 
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(a)  SYMMETRICAL  VORTEX  SHEDDING 
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rj;d  r/*=o 


^;d^/dt^0 


FIG.  2  DIAGRAM  OF  THE  CROSSFLOW  PLANE. 


FIG.  3  DRAG  FORCE  COEFFICIENT  AS  A 
FOR  AN  IMPULSIVELY  STARTED  C 


ASYMMETRIC  MODEL  WITH  SHEDDING 


FIG.  4  DEVELOPMENT  OF  THE  FLOW  FIELD  BEHIND  AN  IMPULSIVELY  STARTED  CYLINDER 
IN  CROSSFLOW. 
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FIG.  6  PITCHING  MOMENT  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK. 
COMPUTED  CURVE  BASED  ON  FREE  PARAMETER  VALUES  USED  FOR  A 
CYLINDER  IN  CROSSFLOW. 
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EXPERIMENTAL  DATA 

O  M=0.5 
A  M=0 .7 


a 

FIG.  7  YAWING  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK.  COMPUTED 

CURVE  BASED  ON  THE  FREE  PARAMETER  VALUES  USED  FOR  A  CYLINDER  IN  CROSSFLOW. 
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EXPERIMENTAL  DATA 


a 


FIG.  8  NORMAL  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK. 
COMPUTED  CURVES  ARE  BASED  ON  A  SET  OF  PARAMETER  VALUES  WHICH 
ARE  MACH  NUMBER  DEPENDENT. 
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EXPERIMENTAL  DATA 

o  <f>  N=o° 

□  *n=90° 

O  $  n=180° 


FIG.  10  YAW  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK  FOR 
M=0.5.  COMPUTED  CURVES  BASED  ON  FREE  PARAMETER  VALUES  WHICH 
ARE  MACH  NUMBER  DEPENDENT. 
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EXPERIMENTAL  DATA 


a 


FIG.  1 1  YAW  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK  FOR 

M-0.7.  COMPUTED  CURVE  BASED  ON  FREE  PARAMETER  VALUES  WHICH  ARE 
MACH  NUMBER  DEPENDENT. 
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EXPERIMENTAL  DATA 


FIG.  12  YAW  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK  FOR 

M=0.9.  COMPUTED  CURVE  BASED  ON  FREE  PARAMETER  VALUES  WHICH  ARE 
MACH  NUMBER  DEPENDENT. 
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EXPERIMENTAL  DATA 


FIG.  13  YAW  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK  FOR  M=1 .1 . 
COMPUTED  CURVE  BASED  ON  FREE  PARAMETER  VALUES  WHICH  ARE  MACH 
NUMBER  DEPENDENT. 
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FIG.  14  COMPARISON  OF  OBSERVED  AND  COMPUTED  VORTEX  PATHS. a=45°  MACHO. 5 


FIG.  15  COMPARISON  OF  OBSERVED  AND  COMPUTED  VORTEX  PATHS 
a  *35°  MACH  =0.5 


- OBSERVED  VORTEX  PATH . 

-  COMPUTED  VORTEX  PATH. 


FIG.  16  COMPARISON  OF  OBSERVED  AND  COMPUTED  VORTEX  PATHS. 
a  =20°  MACH  =  0.5 
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Dtiran  CORRIDOR  IN  WHICH  VORTICES  ARE  OBSERVED. 


FIG.  17  COMPARISON  OF  OBSERVED  AND  COMPUTED  VORTEX  PATHS. 
a  =45°  MACH  =0.9 


FIG.  18  COMPARISON  OF  OBSERVED:  AND  COMPUTED  VORTEX  PATHS. 
a-  35°  MACH  =0.9 


OBSERVED  VORTEX  PATH. 
COMPUTED  VORTEX  PATH. 


FIG.  19  COMPARISON  OF  OBSERVED  AND  COMPUTED  VORTEX  PATHS. 

a  =20°  MACH  =0.9 


NOLTR  73-209 


MODEL 

EXPERIMENTAL  COMPUTED  SOURCE 

L  L.,  R..  DATA  CURVE 

N  N 


12D  4.0  D  0  □  A  PRESENT  STUDY 


Q 


FIG.  20  NORMAL  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK 
FOR  SEVERAL  DIFFERENT  MODELS,  M  =  0.5 
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MODEL 

L  ln 

rn 

EXPERIMENTAL 

DATA 

COMPUTED 

CURVE 

SOURCE 

12D 

4.0  D 

0 

□ 

A 

PRESENT  STUDY 

1  ID 

3.0  D 

0.025D 
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B 

19 

G.  PICK  T 

10D 
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0.250D 
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G.  PICK 
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a 

FIG.  22  YAW  FORCE  COEFFICIENT  AS  A  FUNCTION  OF  ANGLE  OF  ATTACK 
FOR  SEVERAL  DIFFERENT  MODELS,  M  =  0.5. 
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APPENDIX  A 

DERIVATION  OF  PROGRAMMED  EQUATIONS 
The  governing  equations  are: 


ds. 


dP 

1  — - — 

A, 

a* 

il* 

dp 

4.  e  ^2  -  Sox  J 

A* 

d  A  • 

dp 

& 

-S.  j 

In  order  to  define  Xp  c 
augmented  by  the  constraint 
stagnation  points  relative  t 

U>  (  So 

,  )  =  r  dr 

4- 

dr 


growing 

vortices 


r  dP 

Shed  vortices 


(A-l) 


"tcMAct  dli  | 


wtJcx)  r  _JC_  dr  (A-2 ) 

d*  T0  v 


The  dimensionless  velocity  potential,  <t> ,  is  defined  as  follows, 

4>(S)  =  _  <T  -  r*)l  -  tr  V  /jWf 

T  Z-j  L  C  f  -  r  V  f  .  )  - 


•j  =  1 


+  _r_  dr  WCt) 


(A-3) 


Differentiating  Equation  (A-3) ,  the  dimensionless  velocity 
is  determined. 


CON  ly<j#T6  I  -  ?  ( I  +•  )  -  *  r  )  Aj  ("  __J — - 

1  j*  L, 

r  / 


i 


_  ? 

)  (T  -  'VSj)  J 


■f  r  -  dr  J. 
•toM*  d'l  y 


(A-4) 

To  simplify  manipulations,  the  following  functions  are  defined 
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a  (  y,  r j )  --  f  — 1 -  - 

(  ( t  -  5,-) 


(t  -  rVXj-  1 


bit) 


IS 

-  (  i  ♦  f-J  + 

i  e  3 


J 


(A-5a) 


(A-5b) 


Substituting  Equations  (A-5)  into  Equation  (A-4) : 


u>  s  -  i  r  conjugate)  bl3\  +  A.atV,)  +  A  raC^  t8)  +  _J _  dr  | 

*  iom.  d.  di.  ^  * 

(A- 6) 

Substituting  Equation  (A-6)  into  constraints  (A-2) : 

Ai  ft  (  )  t  Xj  a(  J»)  =  -  b  (,  3  ©  t ) 

A  i  ft  +■  A  »  ft  ^  3*-)  ^  ~  b  CT.-v) 

Solving  the  above: 

A  I  &  b  (  So,  )  a,  (  (  3®  I  )  <X  i  "f  <a-aJ  S  %  ) 

ft  ^  f  i )  ft  C  5v)  -  <a-Cf©Sj3,)A(3®i^  3*,) 

As  =  b  (T»,  )«.<?»»,*,)  -  b  ts.-pod?©,,-?,) 

a  C 3o< ^  lx)  ^  C  3| ,)  a  Lf&h  ^ 

Differentiating  the  above  and  substituting  into  the  governing 
equations  produces  differential  equations  for  z .  and  y • .  The 
algebra  involved  is  simplifed  by  defining  the  fallowing  functions. 


(A- 7) 


=  d! 
c 

_  j3  j  (A- 8) 
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T  ( ) 

E  l*,!,*' 
H(3,  1j  ) 


*F  Sj) 
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i.iy-  r V f . ) a  a* 

3  J 
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T/(J  -  rV^j* 
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am  * 

p  i  ?)  _ 

ci 
c  2 

C  3 

C  H  = 
C  S' 

C4, 

C  H 

XII  - 
K  12.  = 

K  13 


~*h 


'  =  3 


(  -L  _  J_  )  +  Y  H  (  f,  1J  )A  J 
d%  $*■  v  *•  ' 


>3=3 
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[  -  *-  C  fo,y  r.)  r  ( s**,  it)  +  *  tfo0f-jr(  to,,  s,)]  /c 
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K  IM 

K  IS’ 
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+  bttox)Ttto,,fa) 


( T,  -  To, ) 
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AX 
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[Ct>  -  (X  t  tQ/  ,s,  )  Cr  C  SoX)  +  T(  t„,Tx)  B(  y*,) 
4  4  /U 


t,  -  t.,) 


K  Ltfj)  =  [  -  atto^fvJCt  t*,y  fj  )  4-ftCtoi,  t,)  Etto^  t3-  )](frfo,)/41 
kM(^)  =  [  -a(  (  f*(j  tt)XFC  Jo?/  ^  )]( T,  ~  1, )  /  4  1 
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k»4  =  [  C7  -  I  bti,|)HU.^ri)  +  p(s4|  )  -b 

-  a. a*,,  h)  Pt5ot  )  ]  /Al  ](  7,-  To,  ) 

K  Z  I  s  C  C  1  ~  b  (.*?»•>)  £  ( 1  )  4  b(?oi)  t  J  (  ^ 

B 1 

KZZ  .-  11  C  2  -  f  ?,)]  (  Si-  *»x  ) 

Bi  31 

K  Z  3  -  C  3  (‘J*-"5-*')  +•  1* 

K  24  r  C  4  (  J,  -  Jo*) 

K  2  S"  =  [  C  S’  -  blJ»x^TC'S*t*  f  *  )  4.  a  (,  ^ot,  5,  )  Q(  ?<tiX  1  ( 

Bl  **  J 

K  2b  r  [  Cb  -  b<  So,)  T  (  ?ox.  y,’)  4.  a^M,  h}  Or  Ulil  ]  (  5»  ~ 

Bi  J 

k  W  (  T  *)  -  [  -  £  c  So*;  f  ^  ^  <xC  f  o»j  S  ^‘ )  t  £  {  T  0)  t  ^'  )  o.  L  Toxj  T ,)  ~](f7~  /  B  j 

KO  ( )  =  [  -  XF(  J.,,  ^  )a  (!•,,  )  +  Xflr,h1j)  OLifo*.,*,)  ](^*-5o,)/Bl 

kZ7  -  t  c  7  -fb  C  -So-^Hft.,,  f,  )  +  <X{  to,,  *,)  PCf.x)  -  b(  fot)  HC*.,,*,  ) 

—  CL  l  J 1#  )  P(S..)]  /  Bl]  <3*-  3.J 

Using  the  above,  explicit  expressions  for  all  derivatives 
can  be  obtained  in  the  following  form: 


ixi 

d% 

„  F  ecx  I  L  uj  ^ 

7 

(A~9a) 

■i  s  3,  ^ -  k 

2^5.  i 
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cl  £ 
d 

^  4  8. 

a7 


V 


dl, 

d*  ' 


[°r'w 


(A-9c) 


The  matrices  D  and  Q  are  defined  as  follows : 


D  (l,  O  =  Re«xl  [  KW  +  K  U] 
DU,Zl  a-  rwa^.tKU  -  kill 

D  <\;3)  s  Reo-V  L  K  13  +  K  1H] 

=  lowa^.lKlM  -K\3] 

D  (  2,  0  =  Rea.1  t  KU+  K»Z] 

D  C?>  1)  =  I  Ml  -  Kl*-] 

D  (?,3  )  =  R«ail  Ki3  tKHl 
0(2,4)  >  lK\3  -KH] 


DC3,0“  Reo-I  LKZl  Hill 
D(3,  2)  =■  L  KU-  Kzi] 

DC3,3)  =  Seal  t  KZ3  +  K2  4] 

DC 3,  4)  =  Jr***].  LkZ4  -K2,3l 

0(4,  l)  =  Re<x\  LkZl  + 

PCH,2)  r  L  KZI-  V<  22-  ] 
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QCi)  » 
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APPENDIX  B 


COMPUTER  PROGRAM  FOR  A  BODY  OF  REVOLUTION  AT 
ANGLE  OF  ATTACK  UP  TO  50  DEGREES 

The  program  listing  and  a  sample  run  are  included  in  this 
section.  The  program  was  developed  on  the  CDC  6400.  Running  time 
is  on  the  order  of  a  minute.  The  sample  run  is  for  a  sharp  nose 
tangent  ogive  body  of  revolution  at  an  angle  of  attack  of 
45  degrees. 


The  program  is  written  to  describe  an  ogive  body  of  revolution 
with  a  sharp  or  blunt  nose.  If  a  sharp  nose  is  to  be  specified, 
variables  listed  by  the  caption  "nose  tip  dimensions"  can  be  ne¬ 
glected.  Allowable  geometries  for  the  blunted  nose  tip  are 
restricted  to  ogive  shapes  of  which  a  circle  is  a  special  case. 


To  use  the  program  the  following  input  is  required. 


Card 


Variables 


Format 


+J 

CO 

P 

e 

CO 

TJ 

U 

cd 

o 

CD 

CO 

CD 

X 

Eh 


U 

o 

m 

Tf  CD 
CD  W 
-P  rd 
rd  o 
0) 

Q^X! 
<D  O 
U  rd 
0) 
CD 


l  6 


N 

A 

XL,  rmax,  ZM2 ,  TL,  OLAP ,  XL1, 

Rn,  ZM1 

XTR,  AKT,  AKL ,  XJf  XST,  S,  Arp/D 
M,  NE,  NGB ,  G,  Y(N+2) 

Ars/D,.GS,  En/D,  AXn,  NNE 


15 

F10. 0 
7F10. 0 

7F10.0 
315,  2F10.0 
4F10.0,  15 


where : 


N 

A 

XL 

rmax 

ZM2 

TL 

OLAP 

XLl 


-  number  of  cases  to  be  run 

-  angle  of  attack  (degrees) 

►  ogive  dimensions 

total  missile  length 


nose  tip  dimensions 


these  parameters  are 
defined  in  the  accom¬ 
panying  figure 


ZM1 
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XTR  fractional  distance  from  missile  nose  to  the  point 

where  transition  from  one  separation  angle  to  another 
takes  place 

AKT  final  separation  angle  (radians) 

AKL  initial  separation  angle  (radians) 

Xj  initial  vortex  strength.  If  not  specified,  the  value 

of  .04  is  assumed. 

XST  nearest  point  to  missile  nose  that  integration  is  allowed 

to  start.  If  not  specified,  the  value  of  .ID  is  assumed. 

S  Strouhal  number 

Arp  initial  radial  perturbation 

M  number  of  integration  steps  between  printing 

NE  error  bound.  Step  size  is  adjusted  to  keep  estimated 

error  between  10~^E  and  io~(NE+3), 

NGB  NGB-1  is  the  number  of  consecutive  steps  to  be  printed 

at  the  start  or  restart  of  integration. 

G  initial  step  size 

Y(N+2)  applies  when  M=0 .  It  is  the  approximate  distance  in 

x  between  printing . 

Ar  first  vortex  is  shed  when  the  radial  location  of  the 

s 

first  two  vortices,  multiplied  by  rmax,  differs  by 
more  than  this  amount . 

rn  this  parameter,  multiplied  by  r,  represents  the  distance 

between  the  flow  separation  point.,  £0-: ,  and  a  newly 
introduced  vortex.  This  distance  is  taken  along  the 
feeding  sheet  of  the  shed  vortex. 

AX  percent  decrease  in  the  circulation  of  a  shed  vortex 

NNE  error  bound  in  effect  after  first  vortex  is  shed 

Values  for  the  parameters  listed  on  this  page  which  are 

used  in  the  current  study  are  found  in  the  right-hand  column  of 

Table  1. 
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The  symbols  on  the  output  page  have  the  following  meanings. 

x  distance  along  body  axis 

R  local  body  radius 

THETAO  0O  ,  0O 

U1  u2 

VORTEX  NUMBER  -  identification  number  for  vortex 


LAM  circulation  strength,  A 

RADIUS 
ANGLE 


cartesian  coordinates  of  vortex 
see  accompanying  figure 
CL  linear  portion  of  normal  force 

CN  nonlinear  portion  of  normal  force 

CML  linear  portion  of  pitching  moment 

CMN  nonlinear  portion  of  pitching  moment 


polar  coordinates  of  vortex 
see  accompanying  figure 


|  CN  =  CL+CN 
|  CM  =  CML+CMN 
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local  normal  force  coefficient, 
yaw  force  coefficient,  Cy 
yawing  moment  coefficient,  CjJj 


CDC 

dx 


£ 


PROGRAM 

VORT  TRACE  CDC  6400  FTN  V3.0-P316  OPT=0 

01/07/74 

14.28.01.  PAGE 

1 

PROGRAM  VORT ( INPUT * OUTPUT *TAPE5=INPUT »TAPE6=0UTPUT) 

VORT 

2 

C. 

..•..MAIN  PROGRAM 

VORT 

3 

COMPLEX  XC 

VORT 

4 

DIMENSION  Y (35) 90(35) 9 XL ( 15) »XC ( 15 ) 9ERROROI) 

VORT 

5 

5 

C0MM0N/GAMM/WRAD*ZMtFLEN»CR*TL»TSL»WRAD2*ZM2tFLEN2>CR2»0LAP 

VORT 

6 

COMMON/BG/N»ORT.PRAD*PRTU» TOP fXLPfPERR»STR. SHED *XAKT*XAKL*XXTR 

VORT 

7 

COMMON/GARB/XL *R»GR »XC * TAL »F VORT 

VORT 

8 

COMMON/ST A/ A *NGB 

VORT 

9 

COMMON/VEL/U 9  IE 

VORT 

10 

10 

COMMON/XLIF/MC9M9CDMN9CDML9XS9CYM9CDCL9CDCN9CY 

VORT 

11 

COMMON/TRAN/XTR  »X  >  AKT  # AKLt ISD 

VORT 

12 

EXTERNAL  DERIV*0UT*TERM2*TERM3 

VORT 

13 

C0NVR=57. 2957795 

VORT 

14 

READ (591005) NDS 

VORT 

15 

15 

DO  21  KT  =  1 9  NDS 

VORT 

16 

READ  (5  9  1 0  0  0 ) A 

VORT 

17 

WRITE (6»200Q ) A 

VORT 

18 

ZM2=0.0 

VORT 

19 

READ (5.1 000 )FLEN*WRAD*ZM*TL*0LAP.FLEN2»WRAD2*ZM2 

VORT 

20 

20 

A=A/CONVR 

VORT 

21 

TAL=TAN (A) 

VORT 

22 

S 

ISEC=0 

VORT 

23 

O 

f 

IF (FLEN2.GT .0.0) ISEC=1 

VORT 

24 

WRITE (692001 ) WRAD9FLEN 9 ZM9TL 

VORT 

25 

25 

CR2=  ( WRAD2*WRAD2*FLEN2*FLEN2) /2./WRAD2 

VORT 

26 

u) 

IF (ISEC.EQ.l) WRITE C6. 2016) WRAD2.ZM2.FLEN2 

VORT 

27 

1 

CR=(WRAD*WRAD+FLEN*FLEN)/2./WRAD 

VORT 

28 

0 

READ (5 9 1000) XTR9AKT9AKL9XLP9XST9STR9PERR 

VORT 

29 

<0 

XXTR=XTR 

VORT 

30 

30 

XAKT=AKT 

VORT 

31 

XAKL-AKL 

VORT 

32 

ISD=0 . 0 

VORT 

33 

IF (AKT.EQ.0.0) AKT=.9251 

VORT 

34 

IF (AKL.EQ.0.0) AKL-13. 12 

VORT 

35 

35 

IF ( AKL.LT .1.6)  ISD=1 

VORT 

36 

IF (XLP.LE.0.0) XLP=.04 

VORT 

37 

IF  (XST.LE. 0 ) XST= .2*WRAD 

VORT 

38 

J=2 

VORT 

39 

IE=1 

VORT 

40 

40 

El=. 000001 

VORT 

41 

E2=E1 

VORT 

42 

DT=. 04 

VORT 

43 

C, 

. DEFINING  COMMON  QUANTITIES 

VORT 

44 

ICOUNT=0 

VORT 

45 

45 

SHED=0  *  0 

VORT 

46 

IF (STR.NE.O) SHED=WRAD/STR/TAL 

VORT 

47 

U=l. 

VORT 

48 

N=5 

VORT 

49 

FVORT=0.0 

VORT 

50 

50 

TOP=0 • 0 

VORT 

51 

DO  10  K=1 9 15 

VORT 

52 

10  XL ( K ) =0 • 0 

VORT 

53 

Ci 

. DETERMINE  INITIAL  INTEGRATION  PARAMETERS 

VORT 

54 

CALL  START (A9DDI9DT9EI9E29XLP9XST) 

VORT 

55 

55 

Ci 

. SET  UP  INTEGRATION 

VORT 

56 

PROGRAM 


VORT 


TRACE 


PAGE 


2 


60 


65 


70 


75 


tu 

f 


80 


85 


90 


95 


100 


105 


110 


VORT  TRACE  CDC  6400 

FTN  V3.0-P316  OPT=0  01/07/74 

14.28.01 

CALL  PROFILE (RtDR*X*DDl?l) 

VORT 

57 

Y(1)=E1*R*C0S(E2> 

VORT 

58 

Y(2)=E1*R*SIN(E2) 

VORT 

59 

Y(3)=-Y(l) 

VORT 

60 

Y ( 4) -Y ( 2) 

VORT 

61 

Y (1)=Y (1) *(1.*PERR> 

VORT 

62 

Y(2)=Y<2)*(1.+PERR) 

VORT 

63 

Y  ( 3)  =  Y ( 3) * ( 1 .-PERR) 

VORT 

64 

Y (4)=Y (4JM1.-PERR) 

VORT 

65 

Y(5)=THETI (DD1) 

VORT 

66 

XS=X 

VORT 

67 

Y ( N+ 1 ) =X 

VORT 

68 

Y (N+3 ) =-l ® 

VORT 

69 

READ(5n0  05)M»NE»NGB9G9Y(N  +  2) 

VORT 

70 

MC=M 

VORT 

71 

M=1 

VORT 

72 

CALL  REST (X9Y9J9G9L9NE9D9O9MC) 

VORT 

73 

CYM=0 • 0 

VORT 

74 

CY  =0 

VORT 

75 

CDCN=SIN(A)*C0S(A)*6.*R*Y(1)*XLP*(1.-R*R/<Y(1)*Y<1)*Y(2)*Y<2) ) ) 

VORT 

76 

1/WRAD/WRAD 

VORT 

77 

CDMN=CDCN*X/2./WRAD 

VORT 

78 

CDCL=2®*SIN (A) *COS (A) *R*R/WRAD/WRAD 

VORT 

79 

XO=CR*CR-FL£N*FLEN 

VORT 

80 

X2=(X0+2.*FLEN*X-X*X)**.5 

VORT 

81 

XO=XO**  » 5 

VORT 

82 

CDML=2.*SIN(A)*C0S(A)*( (WRAD-CR) *< (X+FLEN) *X2/2.-FLEN/2.*X0-CR*CR/ 

VORT 

83 

12.*(ASIN<  (X-FLEN)/CR)-ASIN(-FLEN/CR) ) ) +X*X* (FLEN/2.-X/3. ) )/WRAD/WR 

VORT 

84 

2A0/WRAD 

VORT 

85 

TOP=0 • 0 

VORT 

86 

FVORT=0.0 

VORT 

87 

18  CONTINUE 

VORT 

88 

CALL  FN0L3(U»N*G»L*1»NE*X»Y»D»DERIV»TERM2*0UT»ERR0R) 

VORT 

89 

M-l 

VORT 

90 

CALL  OUT U9Y9D9ERROR9N9L9G) 

VORT 

91 

IF  (X+.OOl.GE.TDGO  TO  20 

VORT 

92 

.....  INTRODUCING  FIRST  NEW  VORTEX 

VORT 

93 

N=7 

VORT 

94 

I  C  =  3 

VORT 

95 

IF (Y (1) **2*Y (2) **2.GT.Y (3) **2*Y (45  **2) GO  TO 

25 

VORT 

96 

WR I TE ( 6  9  2D03 ) 

VORT 

97 

STOP 

VORT 

98 

25  CONTINUE 

VORT 

99 

Y ( 1 0 ) =Y ( 8 ) 

VORT 

100 

Y(9)=Y(7) 

VORT 

101 

Y(8)-Y<6) 

VORT 

102 

Y(7)=Y(5) 

VORT 

103 

Y(6)=Y(2) 

VORT 

104 

Y  <  5 ) =Y ( 1 ) 

VORT 

105 

XL(3)= XL (!)*(!.♦ PR TU) 

VORT 

106 

40  CALL  REST (X9Y9J9G*LfNE9D929MC) 

VORT 

107 

TSL=X+SHED 

VORT 

108 

M=1 

VORT 

109 

CALL  FN0L3 (J.N.G.L.l .NE.X.Y.D.DERIV. TERM3. OUT. ERROR) 

VORT 

110 

M  =  1 

VORT 

111 
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PROGRAM 

VORT  TRACE  CDC  6400  FTN  V3.0-P316  OPT=0  01/07/74 

14.28.01 

CALL  OUT (X*Y»D»ERROR»N*L»H) 

VORT 

112 

IF<X+.001.GE.TL)GO  TO  20 

VORT 

113 

C  .  »  «  s 

•NEW  VORTEX  IS  ADDED 

VORT 

114 

N=N  +  2 

VORT 

115 

115 

IF(N.GT.31)G0  TO  20 

VORT 

116 

Y ( N  +  3) =Y ( N  + 1 ) 

VORT 

117 

Y (N  +  2) =Y (N) 

VORT 

118 

Y ( N+ 1 ) =Y ( N-l ) 

VORT 

119 

Y(N)=Y<N-2) 

VORT 

120 

120 

IC=(N-l)/2 

VORT 

121 

102=101 

VORT 

122 

DO  50  IN=1,IC2 

VORT 

123 

K=2* ( IC-IN) 

VORT 

124 

Y (K+2) =Y (K) 

VORT 

125 

125 

Y (K+l ) =Y (K-l ) 

VORT 

126 

50 

XL(IC+1-IN)=XL(IC-IN) 

VORT 

127 

XL(3)=XL<3)*(1.+PRTU) 

VORT 

128 

I STOP= 1 

VORT 

129 

I COUNT  =  0 

VORT 

130 

130 

GO  TO  40 

VORT 

131 

20 

CNN=CDCL*CDCN 

VORT 

132 

CMM=CDMN+CDML 

VORT 

133 

CSP=CMM/CNN 

VORT 

134 

CYP=0 . 0 

VORT 

135 

135 

IF (DRT.NE.0.0) CYP=CYM/CY 

VORT 

136 

WRITE (6*2004) 

VORT 

137 

WRITE (6.2005) CNN 

VORT 

138 

WRITE(6,2006)CMH 

VORT 

139 

WRITE (6,2007) CSP 

VORT 

140 

140 

WRITE (6,2008) CY 

VORT 

141 

WRITE (6,2009) CYM 

VORT 

142 

WRITE(6,2010)CYP 

VORT 

143 

21 

CONTINUE 

VORT 

144 

STOP 

VORT 

145 

145 

1000 

FORMAT (8F10.0) 

VORT 

146 

1001 

FORMAT (6F 10. 0,15) 

VORT 

147 

1005 

FORMAT (3I5.2F10.0) 

VORT 

148 

2000 

FORMAT ( 1H1 ,40X, 16HANGLE  OF  ATTACK  , F5 . 2 , 2X , 7HOEGREES ) 

VORT 

149 

2001 

FORMAT  (1HO»46X,16HOGIVE  DIMENSIONS/1H  ,6X, 14HMAXIMUM  RADIUS, 6X, 

VORT 

150 

150 

113H  OGIVE  LENGTH, 18X.16HZ  MATCHING  PO I NT , 5X , 1 1 HBOD Y  LENGTH/1H  , 

VORT 

151 

26X,F10.5,14X,F6.3,22X»F10.5,13X,F10.5> 

VORT 

152 

2003 

FORMAT (27HOVORTEX  IDENTITIES  REVERSED) 

VORT 

153 

2004 

FORMAT (1H0»40X,26H*****FINAL  STATISTICS*****) 

VORT 

154 

2005 

FORMAT ( IH  »37X»32HN0RMAL  FORCE  COEFFICIENT=CL*CN=  ,F10.5) 

VORT 

155 

155 

2006 

FORMAT ( IH  »35X»37HPITCHING  MOMENT  COEFFICIENT=CML+CMN=  ,F10.5) 

VORT 

156 

2007 

FORMAT ( 1H  , 41X.23HCENTER  OF  NORMAL  FORCE=  » F 1 0 .5 ) 

VORT 

157 

2008 

FORMAT ( IH  ,42X»22HYAW  FORCE  COEFF ICIENT= » FI  0 .5 ) 

VORT 

158 

2009 

FORMAT ( IH  , 4 1 X , 23HY A W  MOMENT  COEFF I C I ENT= , F 1 0 . 5 ) 

VORT 

159 

2010 

FORMAT  ( IH  ,43X,20HCENTt'R  OF  YAW  F0RCE  =  ,F15.5) 

VORT 

160 

160 

2016 

FORMAT (1HO»46X»15HNOSE  DIMENSI0NS/1H  , 28X  ,  6HR ADIUS , 5X ,  14HMATCHIN 

VORT 

161 

1G  POINT, 5X,6HLENGTH/1H  , 35X , 3 ( F 1 2 . 7 , 3X ) ) 

VORT 

162 

END 

VORT 

163 
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SUBROUTINE 

REST  TRACE  CDC  6400  FTN  V3.0-P316  OPT=0  01/07/74 

14  ®  28 

SUBROUTINE  REST ( X 9 Y , J » 6 9 L 9 NE f D » I T t M ) 

REST 

2 

COMMON/GARB/XL 9R9DR 9 XC9TAL 9 FVORT 

REST 

3 

COMMON/BG/N^DRT^PRADf PRTU 9  TOP 9 XLP 9 PERR 9 S » SHED  9 AKT« AKLsXTR 

REST 

4 

COMMON/GAMM/WRAD 

REST 

5 

5 

DIMENSION  Y(1)«D(1)»XL(15) 9XLR(15) 9XC(15) 

REST 

6 

COMPLEX  XC 

REST 

7 

COMMON/STA/A 

REST 

8 

CONVR=57  *2957795 

REST 

9 

IC= (N-l ) /2 

REST 

10 

10 

IF  ( I T  «  EQ  ®  0 ) GO  TO  10 

REST 

11 

NE=NNE 

REST 

12 

G=XG 

REST 

13 

ANG=Y (N) 

REST 

14 

IF ( ( I C/2) *2®EQ®IC) ANG=3.1415926536-ANG 

REST 

15 

15 

Pl=COS ( ANG) *R 

REST 

16 

P2=SIN ( ANG) *R 

REST 

17 

PAN=AB$( (Y(6)-P2)/(Y(5)-P1) ) 

REST 

18 

PAN=ATAN (PAN) 

REST 

19 

IF (Y (5) «LT  ®P1)PAN=3* 1415926536-PAN 

REST 

20 

20 

Y ( 1 ) =P1 ♦COS (PAN) *PRAD*R 

REST 

21 

Y  (2) =P2+SIN( PAN) *PR AD*R 

REST 

22 

RETURN 

REST 

23 

10  CONTINUE 

REST 

24 

READ (59 1001 ) DRT9XG9PRAD9PRTU9  NNE 

REST 

25 

25 

WRITE (692000 ) XLP9PERR9DRT9PRTU9PRAD9S9SHED9 AKL9AKT9XTR 

REST 

26 

PRTU=PRTU/100. 

REST 

27 

WRITE (692005) 69M9NE9Y (N+2) 9 NNE 

REST 

28 

DRT=DRT*2®*WRAD 

REST 

29 

RETURN 

REST 

30 

30 

1001  FORMAT (4F10®09 15) 

REST 

31 

2000  FORMAT ( 1 HO  9 40X 9 25H*****FREE  PAR AMETERS*****/ 1 H  940X  »24HlNlTI AL  VOR 

REST 

32 

1  TEX  STRENGTH  tF10.5/lH  936X 928HINITIAL  RADIAL  PERTURBATION  9FIO.5/ 

REST 

33 

2 1 H  930X947HRADIAL  ASYMMETRY  AT  WHICH  FIRST  VORTEX  IS  SHED  9FIO.5/ 

REST 

34 

3 1 H  9 28X  9  50HPER  CENT  DECREASE  IN  CIRCULATION  OF  SHED  VORTICES  *F10. 

REST 

35 

35 

45/ 1 H  925X946HL0CATI0N  FACTOR  FOR  NEWLY  INTRODUCED  VORTICES  9FIO.5/ 

REST 

36 

5 1 H  9I8X9I6HSTROUHAL  NUMBER  9 F 1 0 . 5 9 2X 9 32HC0RRESP0ND I NG  SHEDDING  DIS 

REST 

37 

6TANCE  9F10®5/1H  9 1 2X 9 27HSEPARAT I ON  ANGLE  CONSTANTS  9 F 1 0 . 5 9 2X 9 F 1 0  •  5 

REST 

38 

7  9  2X  9  20HPQ I NT  OF  TRANSITION  9F10®5) 

REST 

39 

2005  FORMAT ( 1H0 9 42X 9 1 9H INTEGRAT I ON  OPTIONS/1H  9I8HINITIAL  STEP  SIZE  9FI 

REST 

40 

40 

1  0  ®  7  95X915HPRINT  FREQUENCY  9  15 9 5X 9 1 2HERR0R  BOUND  9 1 5 9 2X 9 1 6HPR I NT  IN 

REST 

41 

2CREMENT  9F10.7/1H  9 42X 9 34HERR0R  BOUND  AFTER  VORTEX  SHEDDING  915) 

REST 

42 

END 

REST 

43 
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FUNCTION 

THETI  TRACE 

CDC  6400  FTN  V3.0-P316  OPT=0  01/07/74 

14.28 

FUNCTION  THETI (DELTA) 

THETI 

2 

COMMON/STA/AL 

THETI 

3 

COMMON/GAMM/WRAD*ZM*FLEN,CR» TL 

THETI 

'  4 

COMMON/TRAN/XTR*X*AKT*AKL* ISD 

THETI 

5 

5 

C0NVR=57.2957795 

THETI 

6 

IF(X/TL.LT.XTR)GO  TO  11 

THETI 

7 

THET I =AKT 

THETI 

8 

RETURN 

THETI 

9 

11 

CONTINUE 

THETI 

10 

10 

IF ( I SD . EQ • 0 ) GO  TO  12 

THETI 

11 

THET I=AKL 

THETI 

12 

RETURN 

THETI 

13 

12 

CONTINUE 

THETI 

14 

THETI=(AKL+C0NVR*DELTA)*(3.13-< 

.116*AL*C0NVR-1 .16) **.5)  THETI 

15 

15 

THETI=THETI/CONVR 

THETI 

16 

RETURN 

THETI 

17 

10 

THETI=3.13*(AKL+C0NVR*DELTA) 

THETI 

18 

THETI=THETI/CONVR 

THETI 

19 

RETURN 

THETI 

20 

20 

END 

THETI 

21 
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SUBROUTINE 


PROFILE  TRACE 


COC  6400  FTN  V3*0-P316  OPT-O  01/07/74  14*28*01 


PAGE 


1 


SUBROUTINE  PROFILE (R^DRsZsDELflC) 

PROFILE 

2 

C0MM0N/GAHM/R02*ZM29XL29CR2fTL9TSL*R0l9ZMl*XLl*CRl90LAP 

PROFILE 

3 

IF ( ICoEGU  1 ) GO  TO  4 

PROFILE 

4 

6 

CONTINUE 

PROFILE 

5 

5 

R0=R02 

PROFILE 

6 

CR=CR2 

PROFILE 

7 

ZM=ZM2«-0LAP 

PROFILE 

8 

IF ( Z*GT • ZH1 ) GO  TO  3 

PROFILE 

9 

RO=RO 1 

PROFILE 

10 

10 

CR=CR1 

PROFILE 

11 

ZM=XL1 

PROFILE 

12 

GO  TO  3 

PROFILE 

13 

4 

SF=TAN ( DEL ) 

PROFILE 

14 

IF  (ZMl®EQoO*0)GO  TO  5 

PROFILE 

15 

15 

Z=XLl-CRl*SF/< <SF*SF+l*)**o5) 

PROFILE 

16 

IF ( ZoLT *ZM1 ) GO  TO  6 

PROFILE 

17 

5 

Z=ZM2+0LAP-CR2*SF/( ( SF*SF+ 1  * ) ** *5) 

PROFILE 

18 

IF  ( Z *LT • ZM1 ) Z-ZM1+  *001 

PROFILE 

19 

GO  TO  6 

PROFILE 

20 

20 

3 

IF (Z»LT «ZM2) GO  TO  20 

PROFILE 

21 

DES=CR*CR- ( ZM-ZM2) **2 

PROFILE 

22 

R=R0-CR+DES***5 

PROFILE 

23 

DR  =  0  «  0 

PROFILE 

24 

IF ( IC • EQ  » 1 ) Z=ZM2 

PROFILE 

25 

25 

RETURN. 

PROFILE 

26 

20 

DES=CR*CR-(ZM-Z) ®*2 

PROFILE 

27 

IF(DES*GEo0*0)GO  TO  30 

PROFILE 

28 

WRITE  (69^0  )ROsZMfZf  DEL  5XLI9CRI 

PROFILE 

29 

STOP 

PROFILE 

30 

30 

30 

DES  =  DES**  ®  5 

PROFILE 

31 

R=RO-CR+DES 

PROFILE 

32 

DR= ( ZM-Z ) /DES 

PROFILE 

33 

RETURN 

PROFILE 

34 

200 

FORMAT (1H  9 25HPR0B ABLE  INPUT  DATA  ERROR  9 6E 1 6 . 8 ) 

PROFILE 

35 

35 

END 

PROFILE 

36 
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SUBROUTINE  START ( AL »DD1 * DT » El * E2 * XLMS * XST ) 

START 

2 

COMMON/GAMM/WRAD«ZM»FLEN»CR*TL 

START 

3 

COMMON/TRAN/XTR»X»AKT  *AKL*ISH 

START 

4 

c 

•  •  •  • 

.INTERNAL  FUNCTIONS 

START 

5 

5 

XR(R1)=((R1**.5-1./(R1**.5>) **2.  ) /4. 

START 

6 

XL (R1»T1*T0)=8.*(XR(R1)+(SIN( <Tl-T0)/2.) ) **2) * <XR <R1 ) 

START 

7 

1 + (COS ( (T1+T0J/2.) >**2)  /C0S(T1)/(R1-1./R1) 

START 

b 

EN1 (D)=(AKL/C0NVR+D)*XK2-ASIN(XK1#TAN(D) ) 

start 

9 

CONVR=57. 2957795 

START 

10 

10 

c 

•  •  •  • 

.DETERMINE  DELTA  CRITICAL 

START 

11 

TAL=TAN ( AL ) 

START 

12 

IF ( ISH. EQ • 0 ) GO  TO  3 

START 

13 

D=ATAN (TAL*SIN (AKL) /1.5) 

START 

14 

GO  TO  10 

START 

15 

15 

3 

CONTINUE 

START 

16 

XK1  =  1 .5/TAN  <AL> 

START 

17 

XK2  =  3 .13 

START 

18 

IF (AL.GT. .17453) XK2=3.13-< . 1 16*AL*C0NVR-1 . 16 )  **.5 

START 

19 

D1=ATAN(1./XK1) 

START 

20 

20 

I F  (  <AKL+D1*C0NVR)*XK2.LT.90.)G0  TO  4 

START 

21 

WRITE (6*2004) 

START 

22 

2004 

FORMAT (21HONO  STARTING  SOLUTION) 

START 

23 

STOP 

START 

24 

4 

CONTINUE 

START 

25 

25 

DDl=Dl/3. 

START 

26 

D=D 1-DD1/2 • 

START 

27 

DP=D 

START 

28 

IC=0 

START 

29 

30 

ER1=EN1 (D) 

START 

30 

30 

IF (ABS(ERl) .LT.E1)G0  TO  10 

START 

31 

IC=IC+1 

START 

32 

IF ( IC.NE. 1 ) GO  TO  20 

START 

33 

ER2=ER1 

START 

34 

DP  =  D 

START 

35 

35 

D=D-DD 1 

START 

36 

GO  TO  30 

start 

37 

20 

DH  =  D 

START 

38 

D=D-ER1* (DP-D) / (ER2-ER1) 

START 

39 

IF (D.GT.D1 ) D= (Dl+QH)/2. 

START 

40 

40 

DP  =  DH 

START 

41 

ER2=ER 1 

START 

42 

GO  TO  30 

START 

43 

10 

CONTINUE 

START 

44 

c 

•  •  •  « 

.TRIAL  AND  ERROR  CALCUATION  FOR  Rl 

START 

45 

45 

CALL  PR0FILE(R*DR*X*D*1) 

START 

46 

IF(X.LT.0.0)X=0.0 

START 

47 

IF  (X/TL.LT.XTR) GO  TO  12 

START 

48 

D=ATAN(TAL*SIN (AKT)/1 .5) 

START 

49 

CALL  PR0FILE(R* DR *X*U*1) 

START 

50 

50 

12 

IF  ( X • GT . XST ) GO  TO  31 

START 

51 

13 

X  =  XST 

START 

52 

CALL  PROFILE ( R *  DR  *  X  *  DEL  *  2 > 

START 

53 

D=ATAN (DR ) 

START 

54 

31 

DDEG=D*CONVR 

START 

55 

55 

D1=D 

START 

56 
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14.28.01 

TG=THETICDI> 

START 

57 

T1=TO+.©01  • 

START 

58 

IT=0 

START 

59 

IP2=1 

START 

60 

60 

TP=T0*1. 047198 

START 

61 

REX=SIN(TO>  -TAM  IIP)  *COS  ( TO) 

START 

62 

TTP=TANCTP) 

START 

63 

70  FE2=FE1 

START 

64 

R1=R£X/«SIIVHT1J-TTP*C0S<T1)  ) 

START 

65 

65 

FE1=ABS  CXL IR1«T1®T0) ) -XLMS 

START 

66 

IF  CABS IFE1 ) .LT.E2)  GO  TO  90 

START 

67 

IT=IT*1 

START 

68 

IF (IT .ME. 1 )  GO  TO  100 

START 

69 

T2=T1 

START 

70 

TO 

T1=T1*DT 

START 

71 

GO  TO  TO 

START 

72 

10©  TH=T1 

START 

73 

T1=T1-F£1* CT2-T1 ) / (FE2-FE1) 

START 

74 

IFCF£1*FE2.LT.0.0) IP2=0 

START 

75 

75 

IF CIP2.E0. 1 ) T1=TH*DT 

START 

76 

J2-JH 

START 

77 

GO  TO  70 

START 

78 

90  E1=R1 

START 

79 

IF CT1. 61. 1. 570796327) Tl=3 . 141592654-T1 

START 

80 

80 

E2=T1 

START 

81 

0D1=01 

START 

82 

RETURN 

START 

83 

END 

START 

84 
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SUBROUTINE 

FN0L3  TRACE.  COC  6400  FTN  V3.0-P316  OPT=0 

01/07/74 

14.28 

SUBROUTINE  FN0L3 ( J 9 NN 9G 9 L 9 MPR 9 NE 9 X 9 Y 9 D 9  DERI  V 9  TERM  9  OUTPUT  9  ERROR ) 

FN0L3 

2 

C 

001 

J=INTEGRATION  METHOD  CONTROL 

FN0L3 

3 

C 

002 

NN=NUMBER  OF  SIMULTANEOUS  DIFFERENTIAL  EQUATIONS 

FN0L3 

4 

C 

003 

G=FIRST  INTERVAL  OF  INTEGRATION 

FN0L3 

5 

5 

C 

004 

L=NUMBER  OF  Y$S  GREATER  THAN  N  TO  BE  PRINTED 

FN0L3 

6 

C 

005 

MPR=PR INT  FREQUENCY 

FN0L3 

7 

C 

006 

NE=CONTROL  FOR  INTERVAL  OF  INTEGRATION 

FN0L3 

8 

C 

007 

X= INDEPENDENT  VARIABLE 

FN0L3 

9 

C 

008 

Y=DEPENDENT  VARIABLE 

FN0L3 

10 

10 

C 

009 

D=ARRAY  CONTAINING  DERIVATIVES 

FN0L3 

11 

C 

010 

DERIV=SUBROUTINE  IN  WHICH  DIFFERENTIAL  EQUATIONS  ARE  FOUND 

FN0L3 

12 

C 

Oil 

TERM=SUBROUTINE  FOR  TERMINATION  CONDITION 

FN0L3 

13 

C 

012 

OUTPUT =SUB ROUTINE  FOR  PRINTING  OUTPUT 

FN0L3 

14 

c 

NRET=TERM  LOOP  COUNTER 

FN0L3 

15 

15 

c 

NPR=COUNT  OF  STEPS  SINCE  LAST  PRINT 

FN0L3 

16 

c 

PC=Y (N+1)=PRINT  CONTROL  OTHER  THAN  STEP  COUNT 

FN0L3 

17 

c 

JJ=J-2=0  FOR  AM 9 - 1  FOR  RK9  +  I  FOR  RK2 

FN0L3 

18 

c 

MK=AM  RK  STEP  COUNT 

FN0L3 

19 

c 

XD=DP  FORM  OF  X 

FN0L3 

20 

20 

c 

DOUBLE  PRECISION  CARD  HAS  BEEN  REMOVED  IN  ORDER  TO  RUN  ON  6400 

FN0L3 

21 

c 

THIS  DECK  RUNS  ON  THE  7090  WITH  OR  WITHOUT  DOUBLE  PRECISION 

FN0L3 

22 

c 

DOUBLE  PRECISION  XD 9  YD 9 YP 9 YC 9 YK *H 9 HC 9 XS * XDS 9 HS »HN * HD24 9 HD6 

FN0L3 

23 

DIMENSION  C  (3) *Y (163) 9  YD ( 163) *YP ( 183) #YC(183> 9D (203) 9 DM (183*4) 

FN0L3 

24 

DIMENSION  DK ( 1 83 9 4 ) ,ERR0R(183> ,YK(183) 

FN0L3 

25 

25 

DATA  EP6*EP11 *M4/l*E-6t l.E-ll*-4/ 

FN0L3 

26 

DATA  (C(K) *K=1*3)/2*,5*I./ 

FN0L3 

27 

DATA  HMAX5/ 1  * E35/ 

FN0L3 

28 

NHTS=0 

FN0L3 

29 

N=NN 

FN0L3 

30 

30 

JJ=J-2 

FN0L3 

31 

H=G 

FN0L3 

32 

HN=H 

FN0L3 

33 

MK=  1 

FN0L3 

34 

NRET=M4 

FN0L3 

35 

35 

EC=Y ( N+3 ) 

FNOL3 

36 

JTEST=1 

FNOL3 

37 

IF  (JJ  «LT a  0)  JTEST=4 

FN0L3 

38 

IF  (NEftEQ.O)GO  TO  15 

FN0L3 

39 

EUP=10ft** (-NE) 

FN0L3 

40 

40 

ELO=EUP* . 00 1 

FN0L3 

41 

EM=EL0*31. 6227766 

FN0L3 

42 

15 

XD=X 

FN0L3 

43 

DO  20  I=I*N 

FN0L3 

44 

ERROR ( I ) =0 • 

FN0L3 

45 

45 

20 

YD ( I ) =Y ( I ) 

FN0L3 

46 

CALL  DERI  V ( X  9  Y  f  D) 

FN0L3 

47 

CALL  TERM  (X9Y9D9F) 

FN0L3 

48 

c 

PRINT 

FN0L3 

49 

50 

CALL  OUTPUT ( X 9 Y »D »ERROR 9 N »L» H) 

FNOL3 

50 

50 

IF  ( NRET )  65960955 

FN0L3 

51 

60 

RETURN 

FNOL3 

52 

55 

WRITE  (693000)  HN 

FN0L3 

53 

3000 

FORMAT (108H1EXECUTION  TERMINATED  BECAUSE  INTERVAL  OF  INTEGRATION 

L  FN0L3 

54 

1ESS  THAN  1.0E  -6  TIMES  INDEPENDENT  VARIABLE  (X).  H  =*1PE15.7) 

FN0L3 

55 

55 

STOP 

FN0L3 

56 
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65 

IF  ( JTEST  •EG*  5  .AND*  H  •  NE® 

HN)  GO  TO  455 

FN0L3 

57 

NPR=0 

FN0L3 

58 

PC=Y (N+l) 

FN0L3 

59 

100 

IF  (JTEST  .EQ.  5  .AND.  H  .NE. 

HN)  GO  TO  455 

FN0L3 

60 

60 

IF  (JJ  .GE.  0)  H=HN 

FNOL3 

61 

IF  (HK  .NE.  0  .OR.  JJ  .NE.  0) 

GO  TO  300 

FN0L3 

62 

c - 

—THE  ADAMS  MOULTON  METHOD 

FN0L3 

63 

200 

HD24=H/24 • 

FN0L3 

64 

JAM  =  0 

FN0L3 

65 

65 

DO  210  1  =  1, N 

FNOL3 

66 

YPI=(55.*DM(I»1) +37.*DM(I.2) >- 

■<59.*DM(I,3> +9.*DM ( I ,4) ) 

FN0L3 

67 

YP(I)=YD(I)+HD24*YPI 

FN0L3 

68 

Y  ( I ) =YP ( I ) 

FNOL3 

69 

210 

CONTINUE 

FN0L3 

70 

70 

X=XD+H 

FN0L3 

71 

CALL  DERIV (X»Y*DM(1»4) ) 

FN0L3 

72 

DO  220  1=1, N 

FNOL3 

73 

YPI=(9.*DM(I»4)*19.*DM(I*1)+DM(I*2) )-5.*DM(I,3) 

FN0L3 

74 

YC(I)=YD(I)*HD24*YPI 

FN0L3 

75 

75 

ERROR (I)=(YP(I)-YC(I) )/14« 

FN0L3 

76 

C  THIS  ADOS  IN  A  2D  CORRECTION 

FNOL3 

77 

YC(I)=YC(I) ♦ERROR ( I ) 

FN0L3 

78 

220 

CONTINUE 

FNOL3 

79 

IF ( NE .NE . 0 ) GO  TO  700 

FN0L3 

80 

td 

l 

80 

GO  TO  455 

FNOL3 

81 

i— 1 

C - 

—THE  RUNGE  KUTTA  METHOD 

FN0L3 

82 

300 

GO  TO  (301,309,308,309,303) , JTEST 

FN0L3 

83 

301 

DO  302  1=1, N 

FN0L3 

84 

YK ( I ) =YD ( I ) 

FN0L3 

85 

85 

302 

CONTINUE 

FN0L3 

86 

XDS=XD 

FN0L3 

87 

GO  TO  309 

FN0L3 

88 

303 

DO  304  1=1, N 

FN0L3 

89 

YK ( I ) =YC ( I ) 

FN0L3 

90 

90 

304 

CONTINUE 

FN0L3 

91 

XDS=XD+H 

FN0L3 

92 

308 

HS=H 

FN0L3 

93 

H=2.*H 

FN0L3 

94 

GO  TO  320 

FN0L3 

95 

95 

309 

X  =  XD 

FN0L3 

96 

JAM=  1 

FN0L3 

97 

DO  310  1=1, N 

FN0L3 

98 

Y  ( I ) =YD ( I ) 

FN0L3 

99 

DK (1,1) =D ( I ) 

FN0L3 

100 

100 

310 

CONTINUE 

FNQL3 

101 

IF  (JTEST  .LE.  2)  CALL  DERIVIX 

»  Y  »  OK ) 

FN0L3 

102 

IF  (MK  .GT.  1  .OR.  JTEST  .GT. 

1)  GO  TO  320 

FN0L3 

103 

DO  315  1=1, N 

FN0L3 

104 

DM ( I ,4) =DK (1*1) 

FN0L3 

105 

105 

315 

CONTINUE 

FN0L3 

106 

320 

DO  335  K=2 , 4 

FN0L3 

107 

HC=H*C ( K-l ) 

FN0L3 

108 

DO  330  1=1, N 

FN0L3 

109 

Y ( I ) =YD ( I )  ♦  HC*DK(I,K-1) 

FN0L3 

110 

no 

330 

CONTINUE 

FN0L3 

111 
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X=XD+HC 

CALL  DERIV(X*Y*DK(1,K> ) 

335  CONTINUE 
HD6=H/6. 

115  DO  340  I  =  1*N 

YPI=DK  (1*1)  +DK(I*4)  +2.*<DK(I,2)  +DK  (1*3)  ) 
YC(I)=YD(I)+HD6*YPI 
340  CONTINUE 

GO  TO  (360*390*370*455*370) *JTEST 
120  360  DO  365  I=1*N 

YP  ( I ) =  YC ( I ) 

365  CONTINUE 
JTEST=3 
GO  TO  308 

125  370  DO  380  1  =  1, N 

YD  ( I ) =YP ( I ) 

YP  ( I ) =YC ( I ) 

380  CONTINUE 
H=HS 

130  XD=XD+H 

JTEST=2 

IF  (MK  .EQ.  1)  GO  TO  309 
GO  TO  451 
390  DO  400  1=1, N 

W  135  ERROR(I)=(YC(I)-YP(I) )/15. 

YC ( I ) =YC ( I ) ♦ERROR ( I ) 

YP ( I ) =YC ( I ) 

400  CONTINUE 
JTEST=5 

140  I F ( NE • NE • 0 ) GO  TO  700 

C - ACCEPT  THE  STEP  SIZE 

450  IF  (JAM  .EQ.  0)  GO  TO  455 
IF  (MK  .EQ.  3  .AND.  JJ  .EQ.  0)  GO  TO  455 
IF  (MK  .NE.  1)  GO  TO  303 

145  IF  ( JTEST  .EQ.  1)  GO  TO  455 

451  DO  452  1=1, N 
Y(I) =YD ( I ) 

452  CONTINUE 
GO  TO  466 

150  455  DO  459  NQ=1*N 

YD ( NQ ) =YC ( NQ ) 

Y ( NQ) =YD ( NQ ) 

459  CONTINUE 

IF  (JJ  .GE.  0)  JTEST= 1 

155  IF  (MK  .NE.  0  .OR.  JJ  .NE.  0  .OR.  NRET  ,NE.  M4)  GO  TO  465 

DO  460  1=1, N 
DM (1*4) =DM (1*2) 

DM (1*2) =DM (1*3) 

DM ( I , 3 ) =DM (1,1) 

160  460  CONTINUE 

465  XD=XD+H 

466  X=XD 

IF  (MK  .EQ.  3)  MK=0 
CALL  DER IV ( X  *  Y  *  D ) 

165  DO  470  I  =  1  * N 
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DM ( I  *  MK  + 1 ) =D { I ) 

FN0L3 

167 

470 

CONTINUE 

FN0L3 

168 

FP=F 

FN0L3 

169 

CALL  TERM  (X*Y*D *F) 

FNQL3 

170 

170 

C - 

—DO  YOU  TERMINATE 

FN0L3 

171 

500 

IF  (ABS(F)  »LE*  EP6 )  GO  TO  800 

FN0L3 

172 

IF  (FP  ®  EQ.  0.)  GO  TO  510 

FNOL3 

173 

IF  (NRET  .NE  •  M4  •OR*  F*FP  «LT.  EP11)  GO 

TO 

805 

FNOL3 

174 

510 

XS  =  XD 

FN0L3 

175 

175 

IF  (MK  .NE.  0  .AND.  H  .EQ.  HN)  MK=MK+1 

FN0L3 

176 

c - 

—DO  YOU  PRINT 

FN0L3 

177 

600 

NPR=NPR ♦ 1 

FN0L3 

178 

IF  (MPR  .GT.  0  • AND »  NPR  .GE.  MPR)  GO  TO 

50 

FN0L3 

179 

IF  (MPR  ,LEo  0  © AND  .  ABS ( Y ( N+ 1 ) -PC )  .GE. 

Y (N+2) ) 

GO  TO  50 

FN0L3 

180 

ISO 

GO  TO  100 

FN0L3 

181 

c - 

--DETERMINING  THE  STEP  SIZE 

FN0L3 

182 

700 

HB  -  HMAX5 

FN0L3 

183 

DO  760  I  =  1 ^ N 

FN0L3 

184 

Z=A8S (ERROR ( I ) ) 

FN0L3 

185 

185 

IF  (YC(I)  .EQ.  0.)  GO  TO  720 

FNOL3 

186 

ZZ  =  ABS(YC(I) ) 

FN0L3 

187 

IF  (EC)  720*710*705 

FN0L3 

188 

705 

IF  (EC  .GT.  ZZ)  ZZ=EC 

FN0L3 

189 

710 

z=z/zz 

FN0L3 

190 

190 

720 

IF (Z.GT.ELO.AND.Z.LT.EUP)  GOTO  750 

FN0L3 

191 

HB  =  AM INI ( HB*EM/ (Z+EP1 1 ) ) 

FN0L3 

192 

GOT0760 

FN0L3 

193 

750 

HB  =  AMIN1 (HB.  1  .  ) 

FN0L3 

194 

760 

CONTINUE 

FN0L3 

195 

195 

IFtHB.EG.l.)  GOTO  450 

FN0L3 

196 

HN=H*HB**.2 

FN0L3 

197 

IF  (MK  .NE.  1)  JTEST=1 

FN0L3 

198 

MK=1 

FN0L3 

199 

IF  (HB.LT.l.)  GOTO  770 

FN0L3 

200 

200 

IF  (ABS(HN)  .GT.  ABS(4.*H))  HN=4.*H 

FN0L3 

201 

NHT  S  =  0 

FN0L3 

202 

GOTO  450 

FN0L3 

203 

770 

HEPS=ABS (X*EP6)  ♦  EP11 

FN0L3 

204 

IF  (ABS(HN)  .LT.  ABS(H/4.))  HN=H/4. 

FN0L3 

205 

205 

IF  (ABS(HN)  .GT.  HEPS)  GO  TO  790 

FN0L3 

206 

NHT S  =  NHTS  ♦  1 

FN0L3 

207 

IF  (NHTS. LT. 10)  GOTO  780 

FN0L3 

208 

JTEST= 1 

FN0L3 

209 

NRET  =  1 

FN0L3 

210 

210 

GOTO  50 

FN0L3 

211 

780 

HN=SIGN (HEPS*HN) 

FN0L3 

212 

790 

IF  (JAM  .EQ.  0)  GO  TO  100 

FN0L3 

213 

DO  795  1  =  1  * N 

FN0L3 

214 

YD ( I ) = YK  (  I ) 

FN0L3 

215 

215 

795 

CONTINUE 

FN0L3 

216 

XD=XDS 

FN0L3 

217 

JTEST=1 

FN0L3 

218 

GO  TO  100 

FN0L3 

219 

C - 

—THE  TERMINATION  LOOP 

FN0L3 

220 

220 

800 

NRET=0 

FN0L3 

221 
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JTEST1=1 

FN0L3 

222 

805 

IF  ( NRET  »NE •  0>  GO  TO  806 

FN0L3 

223 

H=XD-XS 

FN0L3 

224 

GO  TO  50 

FN0L3 

225 

225 

806 

IF(F*FP.LT.O.)  GOTO  810 

FN0L3 

226 

IF  <F*F P2.LT.0.)  GOTO  820 

FN0L3 

227 

GO  TO  800 

FN0L3 

228 

810 

FP2  =FP 

FN0L3 

229 

HP  =H 

FN0L3 

230 

230 

GOTO  830 

FN0L3 

231 

820 

FP  =FP2 

FN0L3 

232 

HP  =H  +  HP 

FN0L3 

233 

830 

NRET=NRET*1 

FN0L3 

234 

H=HP*F/(FP-F) 

FN0L3 

235 

235 

jTEST=4 

FN0L3 

236 

GOTO  300 

FN0L3 

237 

END 

FN0L3 

238 
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SUBROUTINE  TERM2 ( X 9 Y » D * F) 
DIMENSIONY(I) »D(1) 
COMMON/GAMM/WRAD *ZM  ^FLEN^CRsTL 
C0MM0N/BG/N9DRT?PRAD#PRTUf TOP 
5  F=TL-X 

IF (DRT.EQ. 0. 0 ) RETURN 
R1=(Y(1)**2+Y<2)**2)**.5 
R2=<Y(3)**2+Y(4)**2)**.5 
15  IF (ABS<R1-R2) •GTeDRT)F=0«0 
10  RETURN 

END 
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7 

TERM2 
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TERM2 
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TERM2 
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SUBROUTINE  TERM3 ( X * Y * 0 * F ) 

COMMON/BG/N9DRT9PRAD9PRTU9TOP 
COMMON/GAMM/WRAD9ZM9FLEN9CR9TL9TSL 
DIMENSION  YU)»D(1) 

DELT=0.0 

CALL  PROFILE (R*DRf X*DELT «2) 

R2=R*R 

IF<YU)**2+Y<2)**2.GT.R2.AND.Y(3)**2+Y(4)**2.GT.R2)G0  TO  2 
TL=X 

WRITE (692000) 

2  XFIN=TSL 

IF(TL.LT.XFIN)XFIN=TL 

F=XFIN-X 

RETURN 

2000  FORMAT (51H0EXECUT ION  TERMINATED  DUE  TO  NUMERICAL  INSTABILITY.) 
END 


TERM3 

2 

TERM3 

3 

TERM3 

4 

TERM3 

5 

TERM3 

6 

TERM3 

7 

TERM3 

8 

TERM3 

9 

TERM3 

10 

TERM3 

11 

TERM3 

12 

TERM3 

13 

TERM3 

14 

TERMS 

15 

TERMS 

16 

TERM3 

17 
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SUBROUTINE 

DERI V  TRACE  CDC  6400  FTN  V3.0-P316  OPT  =  0 

01/07/74 

14.28 

SUBROUTINE  DERIV(X*Y*D> 

DERI  V 

2 

COMPLEX A, B »T»E *XF » G» XZ.XF 1  * XF2 , XF3 . DZOP .X J *  X JC *  Cl • C2 > C3»C4. C5 * 

DERI  V 

3 

1C6*K11»K12*KI3»K14*K21*K22*K23*K24*K15»K16*K25*K26»DTH»0Z0*TH2* 

DERIV 

4 

2KN1»KN2»WN»XK.XKK*XWN*CIR*XC*L1.L2>X1*X2*X0*X0P»P*H»C7»K19.K29» 

XNP  DERIV 

5 

5 

INTEGER  ZO9ZOP9ZI ?Z2®ZJ 

DERIV 

6 

COMMON/GARB/XLR*R.DR.XC*TAL*FVORT 

DERIV 

7 

COMMON/ BG/N  9DRT9PRAD9PRTU9TOP 

DERIV 

8 

COMMON/VEL/U 9  IE 

DERIV 

9 

DIMENSION  Y ( 1 ) »D(1) *WN(15) »DD(4»4) tQ(4,l) »XL<15) *XC(15) ,XZ<2) * 

DERIV 

10 

10 

1A(2.2).B(2).T(2»15)»E(2»15)»XF<2.15>,G(2)»P<2)»H(2»2)  »XNPU5> 9 

DERIV 

11 

3XLR ( 15 ) 

DERIV 

12 

CALL  PROFILE ( R t DR  $ X 9 DELTA  *  2) 

DERIV 

13 

XK=CMPLX(0«OfI®0) 

DERIV 

14 

XKK=XK*R 

DERIV 

15 

15 

IC=(N-l)/2 

DERIV 

16 

DO  22  KQ= 1 9 I C 

DERIV 

17 

22  XL(KO)=XLR(KQ)/R 

OERIV 

18 

C. 

..OBTAIN  THE  STAGNATION  POINTS 

DERIV 

19 

THETA=Y(N)*(l.-TOP> 

DERIV 

20 

20 

Sl=R*COS( THETA) 

DERIV 

21 

S2=R*SIN(THETA) 

DERIV 

22 

XZ(1>=  CMPLX (SI »S2) 

DERIV 

23 

S1=-S1 

DERIV 

24 

XZ (2) =CMPLX (SI9S2) 

DERIV 

25 

25 

IF<IC.EQ.2.0R.IC.NE.2*<IC/2) )GO  TO  6 

DERIV 

26 

DTH=XZ (1) 

DERIV 

27 

XZ ( X ) =XZ ( 2 ) 

DERIV 

28 

XZ ( 2) =  DTH 

DERIV 

29 

6  CONTINUE 

DERIV 

30 

30 

C. 

...•TRANSFORMS  VORTEX  POSITIONS  TO  COMPLEX  COORDINATES 

DERIV 

31 

DO  5  1=1*10 

DERIV 

32 

5  XC( I)=CMPLX<Y(2*I-1) *Y(2*I) > 

DERIV 

33 

c. 

...•CALCULATE  FUNCTIONAL  VALUES 

DERIV 

34 

DO  4  L  =  1 9  2 

DERIV 

35 

35 

XI =XZ ( L ) 

DERIV 

36 

B (L) = (R*R/Xl/Xl+1 ® )/R+FVORT/Xl 

DERIV 

37 

G (L)=-2.*R/X1/X1/X1-FV0RT/X1/X1 

DERIV 

38 

P (L)=DR*  < 1./X1/X1-1./R/R) 

DERIV 

39 

DO  4  K= 1 9  I C 

DERIV 

40 

40 

X  J  =  XC ( K ) 

DERIV 

41 

XJC=CONJG(XJ) 

DERIV 

42 

XF1=1./ (X1-XJ)-1./(X1-R*R/XJC) 

DERIV 

43 

XF2=1./(X1-R*R/XJC) **2-1./ (Xl-XJ) **2 

DERIV 

44 

XF3=-DR*2.*R/< (X1-R*R/XJC)**2) /XJC 

DERIV 

45 

45 

XF (L*K)=<R/ <X1-R*R/XJC) /XJC) **2 

DERIV 

46 

E(L9K)=1./(X1-XJ)**2 

DERIV 

47 

IF (K.GT.2) GO  TO  3 

DERIV 

48 

A(L*K)=XF1 

DERIV 

49 

T ( L  9  K ) =  XF 2 

DERIV 

50 

50 

H  ( L  9  K ) =XF3 

DERIV 

51 

GO  TO  4 

DERIV 

52 

3  B(L)=B(L)+XL(K)*XF1 

DERIV 

53 

G  CL) =  G  ( L ) +  XL (K) *XF2 

DERIV 

54 

P(L)=P(L)+XL(K)*XF3 

DERIV 

55 

55 

4  CONTINUE 

DERIV 

56 
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65 


70 


75 


85 


90 


95 


100 


105 


110 


C . CALCULATE  THE  TWO  UNKNOWN  CIRCULATIONS 

ZO=l 
ZOP  =  2 
Z  1  =  1 
Z2  =  2 

Al=REAL(-B(Z0)*A(Z0P*Z2) ♦B(Z0P)*A(Z0*Z2) ) 

B1=REAL <-B (ZOP) *A (Z0*Z1) +  B <ZO)*A <Z0P»Z1 ) ) 

C  =  REAL (  A ( ZO  *  Z1 ) *A (  ZOP  *  Z2 ) -A ( ZO  *  Z2 ) *A ( ZOP  *  Z 1 ) ) 

XL(l)=Al/C 

XL (2) =B1/C 

XLR(1)=XL(1)*R 

XLR ( 2 ) =XL ( 2 ) *R 

IF (U.LT.O) RETURN 

C. ....CALCULATE  THE  FLUID  VELOCITY  AT  EACH  VORTEX  LOCATION 
DO  20  K=1,IC 
Xl=XC (K ) 

XWN=R*DR/X1/TAL-XK*(1.+R*R/X1/X1)-XK*FV0RT/X1 
DO  10  L=mc 
X J=XC ( L ) 

IF ( L .NE .K ) GO  TO  8 

XNP(L) =XKK*XL (L)/ < X 1 -R*R/CON JG ( X J ) ) 

GO  TO  10 

8  XNP(L)=-XKK*XL (L) * ( 1 ./ (Xl-XJ) -1 ./ ( Xl-R*R/CONJG (XJ) ) ) 

10  CONTINUE 

WN (K) =XWN+XNP (K) 

DO  15  MBS=1*IC 
I F ( MBS . EQ • K ) GO  TO  15 
WN ( K ) =WN ( K ) +  XNP ( MBS ) 

15  CONTINUE 

20  WN ( K ) =CON JG ( WN ( K ) ) 

C . DETERMINE  THE  DERIVATIVE  OF  THETA 

DELT  A  =  ATAN ( DR ) 

THETA  I  =  THET I (DELTA) 

D(N)=(.7*TAL/R)*(THETAI-Y(N) ) 

C . CALCULATE  DERIVATIVES  FOR  ALL  FREE  VORTICES 

IF ( IC.LT.3) GO  TO  40 
DO  30  IL=3 »IC 
RWN=REAL (WN ( IL) ) 

XIWN=AIMAG (WN ( IL) ) 

D(2*IL)=XIWN*TAL 
30  D(2*IL-1)=RWN*TAL 
40  CONTINUE 

C . CALCULATE  DD  MATRIX  FOR  ATTACHED  VORTICES 

X2=XC (2) 

X1=XC(1) 

XO=XZ  (1) 

X0P=XZ<2) 

C 1 =- ( A (ZOP*Z2>*E(Z09Z1)-A(ZO*Z2)*E<ZOP*Z1> )/C 
C2=-(A (ZOP«Z2)*XF (ZOtZl)-A (Z0,Z2)*XF (Z0P*Z1) ) /C 
C3=-(A(Z09Z1)*E(Z0P?Z2)-A(Z0P»Z1)*E(Z0,Z2> ) /C 
C4=-(A (Z09Z1)*XF(Z0P»Z2)-A(Z0P,Z1)*XF(Z0*Z2) ) /C 
K11=C1* (Xl-XO ) ♦ 1 • 

K12=C2* (Xl-XO) 

Kl3=(C3-(B(ZO)*E(ZOp9Z2)-B(ZOP)*E(ZO»Z2) ) /A1 ) * (Xl-XO) 
K14=(C4-(B(Z0)*XF (ZOP*Z2)-8 (ZOP) *XF(Z0?Z2) ) /Al ) * (Xl-XO) 
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76 
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81 
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83 
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84 
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85 
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86 

DERIV 

87 

DERIV 

88 

DERIV 

89 

DERIV 

90 

DERIV 

91 

DERIV 

92 

DERIV 

93 

DERIV 

94 

DERIV 

95 

DERIV 

96 

DERIV 

97 

DERIV 

98 

DERIV 

99 

DERIV 

100 

DERIV 

101 

DERIV 

102 

DERIV 
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DERIV 
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DERIV 
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DERIV 
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115 


120 


125 


130 


135 


140 


145 


150 


155 


160 


165 


K21=(C1-(B(Z0P)*E(Z09Z1)-B(Z0)*E(Z0P,Z1> ) /B 1 ) * ( X2-XOP )  DERIV 

K22=(C2-(B(ZOP)*XF(Z09Z1)“B(ZO)*XF(ZOP,Z1) ) /B 1 ) * ( X2-X0P )  DERIV 

K23=C3* (X2-X0P) +1 .  DERIV 

K24=C4* (X2-XOP)  DERIV 

DD (1 91)=REAL (K11+K12)  DERIV 

DD  (1  92)  =AIHAG  (K12-KH)  DERIV 

DD (1*3)=R£AL(K13+Kl4>  DERIV 

DD (1*4)=AIMAG(K14-K13)  DERIV 

DD (29 1 )=AIMAG(K11+K12)  DERIV 

DD ( 2  9  2 ) =REAL (K11-K12)  DERIV 

DD (2*3)=AIMAG(K13+K14)  DERIV 

DD (2*4)=REAL (K13-K14)  DERIV 

DD (391)=REAL(K21+K22)  DERIV 

DD (392)=AIMAG(K22-K21)  DERIV 

DD (3?3)=REAL ( K23+K24 )  DERIV 

DD (394) =AIHAG (K24-K23)  DERIV 

DD  <4*1) =AIMAG (K21+K22)  DERIV 

DD (4*2) =REAL (K21-K22)  DERIV 

DD (4*3) =AIMAG(K23*K24)  DERIV 

DD (4*4) =REAL (K23-K24)  DERIV 

C..... CALCULATE  THE  Q  VECTORS  DERIV 

C5=  <-A (Z0P*Z2) *T (ZO,Zl) ♦ A ( ZOP * Z 1 ) *T ( ZO * Z2 )  ) /C  DERIV 

C6  =  (  A(ZO*Zl)*T  (Z0P,.Z2>-A(Z09Z2)*T<Z0P9Z1)  ) /C  DERIV 

C7=(A(Z0P9Z1)«H(Z09Z2) + A ( ZG 9 Z2 ) *H ( ZOP * Z 1 ) -A ( ZOP 9 Z2 ) *H ( ZO 9 Z 1 >  DERIV 

1-A ( ZO  *  Z1 ) *H ( ZOP  t  Z2) )/C  DERIV 

K 15=  (C5-(A  (Z0P9Z2) *G(ZO)-B (ZOP)*T (Z0*Z2)  )/Al)*(Xl-Xt>)  DERIV 

K16= (C6*  (  T  <Z0P9Z2) *B(ZO) -A (Z09Z2) *G (ZOP) ) /A  1 ) * ( X 1 -XO )  DERIV 

K19=(C7* ( A ( ZO  9  Z2 ) *P (ZOP ) ♦B ( ZOP ) *H ( ZO * Z2 ) -A ( ZOP » Z2 ) *P ( ZO ) -B ( ZO ) *  DERIV 
1H ( ZOP  9  Z2 ) )/Al)*<Xl-XO)  DERIV 

K25=  ( C5+ (  A(ZOPtZl)*G(Z0)-T(Z09Zl)*B(Z0P) ) /B 1 ) * { X2-X0P )  DERIV 

K26=(C6*(  A (ZO9ZI) *G (ZOP) -8 (ZO) *T  <Z0P9Z1) ) /B1 ) * ( X2-X0P )  DERIV 

K29=(C7«MA(Z0P»Z1)*P(Z0) +& ( ZO ) *H ( ZOP 9 Z 1 ) -A ( ZO 9 Z 1 ) *P ( ZOP ) -B ( ZOP ) *  DERIV 

lH(ZOtZl) ) /B 1 ) * ( X2-XOP )  DERIV 

DT2=D (N)*R* (l.-TOP)  DERIV 

DTH  =  CMPLX  <DR*DT2)  DERIV 

TH2  =  CMPLX (0.0 9  THETA)  DERIV 

DZO=CEXP (TH2) *DTH  DERIV 

DTH=CMPLX (-DR9DT2)  DERIV 

DZ0P=-CEXP (-TH2) *DTH  DERIV 

IF(IC.EQ.2.0R.IC.NE.2*(IC/2) ) GO  TO  55  DERIV 

TH2=DZO  DERIV 

DZO=DZOP  DERIV 

DZ0P=TH2  DERIV 

55  CONTINUE  DERIV 

L1=-K15*DZ0-K16*DZ0P*TAL*WN(1)-DR* < Xl-XO) /R-K19  DERIV 

L2  =  -K25*DZ0-K26*DZ0P  +  TAL*WM2)  -DR*  (X2-XOP)  /R-K29  DERIV 

IF ( IC.LT.3) GO  TO  50  DERIV 

DO  60  M  =  3  9  I C  DERIV 

DTH=CMPLX (D <2*M-1 ) 9D(2*M) )  DERIV 

TH2=CONJG(DTH)  DERIV 

ZJ=M  DERIV 

C I R  =  XL ( M )  DERIV 

KN1= (Xl-XO) *  <-A (Z0P,Z2) *E (ZO9ZJ) ♦ A ( ZO 9 Z2 ) *E ( ZOP 9 Z J ) ) /A1*CIR  DERIV 

KN2= (XI -XO ) * (-A ( ZQP9Z2) *XF ( ZO  9  Z J ) ♦ A ( ZO 9 Z2 ) *XF ( ZOP 9 Z J )  ) /A1*CIR  DERIV 

L1=L1-KN1*DTH-KN2*TH2  DERIV 
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SUBROUTINE 

DERI V  TRACE  CDC  6400  FTN  V3.0-P316  OPT  =  0 

01/07/74 

14.28 

KN1=(-E<Z0P*ZJ>*A<Z09Z1) +A (Z0P»Z1)*E (ZO9ZJ) )*<X2-X0P)*CIR/61 

DERI  V 

167 

KN2=(-XF(ZOP9ZJ)*A(ZO  9  Z1 ) +XF (ZO>ZJ)*A(ZOP»Zl) ) * ( X2-X0P ) I R/B 1 

DERI  V 

168 

60 

L2=L2-DTH*KN1-TH2*KN2 

DERI  V 

169 

50 

Q  ( 1  9  1  >  =  REAL (LI ) 

DERIV 

170 

170 

Q(2»1>=AIMAG(L1) 

DERI  V 

171 

Q  (3  v 1 ) =  REAL ( L2) 

DERIV 

172 

q{4,1)=AIMAG(L2) 

DERIV 

173 

11=4 

DERIV 

174 

IX=l 

DERIV 

175 

175 

CALL  MATRIX ( DD *  1 1  *  1 1  * Q *  IX  *  DETERM  *  ID) 

DERIV 

176 

D  ( 1 ) =Q ( 1  *  1 ) 

DERIV 

177 

D  ( 2 ) =0 ( 2  *  1 ) 

DERIV 

178 

D  ( 3 ) =Q ( 3  *  1 ) 

DERIV 

179 

D  ( 4 ) =Q ( 4  » 1 ) 

DERIV 

180 

180 

RETURN 

DERIV 

181 

END 

DERIV 

182 

w 
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i 
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25 


30 


35 
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45 


50 


55 


SUBROUTINE  OUT ( X  * Y » D t ERROR t N #L * H ) 

COMPLEX  XC » F 1 

DIMENSION  XC ( 15) *  XL ( 15) ♦ D ( 1 )  ,  ERROR ( 1 ) » Y (1) *XLR (15) 
COMMON/ ST  A/ A  ,  NGB 

COMMON/XL IF/MC,M , CDMN t CDML , XS , CYM , CDCL1 *CDCN1 ,CY1 

COMMON/GARB/XLR»R,DR,XC,TAL 

COMMON/VEL/U *  IE 

COMMQN/GAMM/RO  ,  ZM  ,  FLEN  , CR  »  TL 

COMMON/BG/IBLANK  »DRT 

C0NVR=57. 2957795 

I F ( DRT • NE • 0 • 0  ) GO  TO  6 

Y(1)=(Y(1)-Y(3) )/2. 

Y(3)=-Y(l> 

Y(2)=(Y<2)*Y(4>)/2. 

Y (4) =Y (2) 

8  CONTINUE 
U=-U 

CALL  DERI V (X*  Y »D) 

U=-U 

C . CALCULATE  THE  FORCE  COEFFICIENTS 

CDCL=2.*SIN (A) *COS (A) *R*R/RO/RO 
F1  =  CMPLX  (0*0,0. 0) 

IC=(N-l)/2 
DO  22  KQ= 1 , I C 
22  XL(KQ)=XLR(KQ)/R 
DO  10  1  =  1, IC 

10  F1=XL ( I ) * (XC ( I)-R*R/CONJG(XC (I) ) )+Fl 
Fl=4.*SIN(A)*COS(A) *Fl*R/RO/RO 
CDCN=REAL (FI) 

CY  =  A IMAG  (FI ) 

CD=0 • 0 

IF (X.NE.XS) CD= (CDCN-CDCN1) / (X-XS) 
CD=CD*3.1416*R/2./SIN(A)/SlN(A) 

C . CALCULATE  MOMENT  COEFFICIENTS 

CDML= (X+XS)*(CDCL-CDCLl)/4./RO  +CDML 
CDMN=(XS+X) * (CDCN-CDCN1 )/4*/RO  +  CDMN 

CYM= (X+XS ) * ( CY-CY 1 ) /4./R0  +  CYM 

C . REINITIALIZE  FORCE  COEFFICIENTS 

CDCL 1 =CDCL 
CDCN 1 =CDCN 
C Y 1 =CY 

IF (X.EQ.XS) CIT=0 
XS  =  X 
M=M- 1 

IF (X.EQ.TL) GO  TO  16 
CIT=CIT+1 

IF (CIT.LT.NGB) GO  TO  16 
I F ( MC • EQ • 0 ) GO  TO  15 
IF ( M • EQ • 0 ) GO  TO  16 
RETURN 

15  IF ( M • EQ • 0 ) GO  TO  17 

IF (X-Y (N+l) *LT . Y (N+2) ) RETURN 
Y (N  +  l ) =X 
GO  TO  17 

16  M=MC 
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OUT 

TRACE  CDC  6400  FTN  V3.0-P316  OPT=0  01/07/74 

14.28.01 

17 

DDEG=CONVR*Y (N) 

OUT 

57 

HPP=0 • 0 

OUT 

58 

IF ( ABS { H) • LT • 1 • ) HPP=H 

OUT 

59 

WRITE (692000) X9R9DDEG9HPP 

OUT 

60 

60 

WRITE (6*2004) 

OUT 

61 

DO  20  KT=1 9 IC 

OUT 

62 

K  =  1 +  IC-KT 

OUT 

63 

IF (IC.EQ.2)K=KT 

OUT 

64 

XR2= (Y (2*K-1)**2+Y (2*K> **2)**.B 

OUT 

65 

65 

ANG2=C0NVR*ATAN ( Y ( 2*K-1 ) /Y < 2*K) ) 

OUT 

66 

I F  (  Y  ( 2*K ) .LT.0.0) ANG2=Y(2*K-1)/ABS(Y(2*K-1) ) * ( ABS ( ANG2) +90.0) 

OUT 

67 

20 

WRITE (6 9 2001) KT 9  XL (K) 9 XR2 9 ANG2 9 Y ( 2*K-1 ) 9Y (2*K) 

OUT 

68 

WRITE(692002)CDCL9CDML9CDCN9CDMN9CD9CY9CYM 

OUT 

69 

IF(I£.EQ.O)RETURN 

OUT 

70 

70 

WRITE (6,2003) (ERROR (II) 9 I 1=1 9 N) 

OUT 

71 

RETURN 

OUT 

72 

2000 

FORMAT (5H0X  =  9 F 1 2 . 7 9 4X 9 4HR  =  9 F 1 2 . 7 9 4X 9 9HTHETA0  =  9 F 1 0  * 5 9 4X 9 1 2HST 

OUT 

73 

1EP  SIZE  *  9E14.Q) 

OUT 

74 

2001 

FORMAT (1H  920X9 1396X95 (F14.793X) ) 

OUT 

75 

75 

2002 

FORMAT ( 1H  96X92HCL9I2X94H  CML9  11X9  2HCN9  12X93HCMN9  HX9  2HCD9l2Xt2HCY 

OUT 

76 

I912X93HCYM9/7F14.7) 

OUT 

77 

2003 

FORMAT (  8H  ERRORS  98E14.6) 

OUT 

78 

2004 

FORMAT (1H  9l5X914HVORTEX  NUMBER  9 7X 9 3HLAM 9 13X 9 6HRADIUS* 1 IX 9 5HANGLE 

OUT 

79 

I9I3X9IHY9I6X9IHZ) 

OUT 

80 

80 

END 

OUT 

81 
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SUBROUTINE 

MATRIX  TRACE  CDC  6400 

FTN  V3 «  0-P316  OPT  =  0 

01/07/74 

14«28»01«  PAGE 

SUBROUTINE  MATRIX < A , N1 *  1 1 > B *  Ml , DETERM » ID ) 

MATRIX 

2 

C 

MATRIX  INVERSION  WITH  ACCOMPANYING  SOLUTION 

OF  LINEAR  EQUATIONS 

MATRIX 

3 

C 

MATRIX 

4 

C 

TEST  FOR  LOST  OF  DIGITS  DUE  TO  SUBTRACTION 

MATRIX 

5 

5 

C 

MATRIX 

6 

DIMENSION  A  ( 1 1 9 1 1 ) 'B(Il'Ml) 9 INDEX (5093) 

MATRIX 

7 

EQUIVALENCE  (IROW9 JROW)  *  (I COLUM  t  JCOLUM  )  »  ( AM  AX  9  T9  SWAP) 

MATRIX 

8 

DATA  EPS4/1.E-7/ 

MATRIX 

9 

c 

MATRIX 

10 

10 

ID=1 

MATRIX 

11 

M=M1 

MATRIX 

12 

N=N1 

MATRIX 

13 

10 

DETERM= 1 «  0 

MATRIX 

14 

15 

DO  20  J=1 

MATRIX 

15 

15 

20 

INDEX (J93)  =  0 

MATRIX 

16 

30 

DO  550  1=1 *N 

MATRIX 

17 

c 

MATRIX 

18 

c 

SEARCH  FOR  PIVOT  ELEMENT 

MATRIX 

19 

c 

MATRIX 

20 

20 

40 

AMAX=0 • 0 

MATRIX 

21 

45 

DO  105  J=1*N 

MATRIX 

22 

IF  ( I NDEX ( J 9 3 )  ,EQ8  1)  GO  TO  105 

MATRIX 

23 

60 

DO  100  K=1 9  N 

MATRIX 

24 

IF ( INDEX ( K  *3 ) -1 )  80  9  100?  100 

MATRIX 

25 

25 

80 

IF  (AMAX  • GE  ®  ABS ( A ( J  9  K )  )  )  GO  TO  100 

MATRIX 

26 

85 

I ROW  =  J 

MATRIX 

27 

90 

I COLUM=K 

MATRIX 

28 

AMAX  =  ABS (A ( JtK) ) 

MATRIX 

29 

100 

CONTINUE 

MATRIX 

30 

30 

105 

CONTINUE 

MATRIX 

31 

INDEX ( I COLUM 9 3 )  =  INDEX  < ICOLUMf 3)  +1 

MATRIX 

32 

260 

I NDEX ( I 9 1 ) = I ROW 

MATRIX 

33 

270 

INDEX(I92)=IC0LUM 

MATRIX 

34 

IF (BMAXaEQaAMAX) ID-2 

MATRIX 

35 

35 

c 

MATRIX 

36 

c 

INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

MATRIX 

37 

c 

MATRIX 

38 

130 

IF  (IROW  » EQ 0  ICOLUM)  GO  TO  310 

MATRIX 

39 

140 

DETERM=~DETERM 

MATRIX 

40 

40 

150 

DO  200  L  =  1 9  N 

MATRIX 

41 

160 

SWAP=A (IROWfL) 

MATRIX 

42 

170 

A(IROW9L)=A(ICOLUM9L) 

MATRIX 

43 

200 

A<ICOLUMfL)=SWAP 

MATRIX 

44 

IF (M)  3109  3109  210 

MATRIX 

45 

45 

210 

DO  250  L=l9  M 

MATRIX 

46 

220 

SW AP=B ( IROW9L) 

MATRIX 

47 

230 

B  < I RO W  9  L ) =B ( I COLUM  9  L ) 

MATRIX 

48 

250 

B<ICOLUMfL)=SWAP 

MATRIX 

49 

c 

MATRIX 

50 

50 

c 

DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

MATRIX 

51 

c 

MATRIX 

52 

310 

PIVOT  =A (ICOLUMt ICOLUM) 

MATRIX 

53 

DETERM=DETERM*PIVOT 

MATRIX 

54 

330 

A (ICOLUMfICOLUM)=l*0 

MATRIX 

55 

55 

340 

DO  352  L  =  1 9  N 

MATRIX 

56 
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SUBROUTINE 

MATRIX  TRACE.  CDC  6400  FTN  V3.0-P316  OPT=0 

01/07/74 

14.28.01 

IF  (PIVOT  *NE.  0.  )  GO  TO  350 

MATRIX 

57 

A  ( ICOLUM , L )  =  0. 

MATRIX 

58 

GO  TO  352 

MATRIX 

59 

350 

A  ( ICOLUM  *  L ) =A ( ICOLUMt L> /PIVOT 

MATRIX 

60 

60 

352 

CONTINUE 

MATRIX 

61 

355 

IF ( M )  380  *  380  ,  360 

MATRIX 

62 

360 

DO  375  L=1,M 

MATRIX 

63 

IF  (PIVOT  .NE*  0*  )  GO  TO  370 

MATRIX 

64 

B ( ICOLUM  , L )  =  0. 

MATRIX 

65 

65 

GO  TO  375 

MATRIX 

66 

370 

B  ( ICOLUM 9 L) =B ( ICOLUM? L) /PIVOT 

MATRIX 

67 

375 

CONTINUE 

MATRIX 

68 

C 

MATRIX 

69 

C 

REDUCE  NON-PIVOT  ROWS 

MATRIX 

70 

70 

C 

MATRIX 

71 

BMAX=0  «  0 

MATRIX 

72 

380 

DO  550  L1=1?N 

MATRIX 

73 

390 

IF  <L1  .EQ.  ICOLUM)  GO  TO  550 

MATRIX 

74 

400 

T=A (LI  , ICOLUM) 

MATRIX 

75 

75 

420 

A (LI ? ICOLUM) =0. 0 

MATRIX 

76 

430 

DO  450  L= 1 ?  N 

MATRIX 

77 

SUB  =A ( ICOLUM  *  L )  *T 

MATRIX 

78 

A ( L 1 »  L ) =A (LI , L ) -SUB 

MATRIX 

79 

C 

MATRIX 

80 

60 

C 

MATRIX 

81 

IF ( INDEX ( L 1  * 3 ) .EQ • 1 )  GO  TO  450 

MATRIX 

82 

I F ( ABS ( A ( L 1 ?  L ) )-EPS4*  ABS(SUB) ) 449, 449? 450 

MATRIX 

83 

449 

BMAX=AMAX1 (BMAX,ABS(A(L1,L) ) ) 

MATRIX 

84 

450 

CONTINUE 

MATRIX 

85 

85 

455 

IF  (M )  550,  550,  460 

MATRIX 

66 

460 

DO  500  L=1,M 

MATRIX 

87 

500 

B (LI ,L) =B (LI ,L) -B  < ICOLUM, L) *T 

MATRIX 

88 

550 

CONTINUE 

MATRIX 

89 

C 

MATRIX 

90 

90 

C 

INTERCHANGE  COLUMNS 

MATRIX 

91 

C 

MATRIX 

92 

600 

DO  710  1=1, N 

MATRIX 

93 

610 

L=N+ 1 -I 

MATRIX 

94 

620 

IF  (INDEX (L,l)  « EQ »  INDEX(L*2))  GO  TO  710 

MATRIX 

95 

95 

630 

JROW=INDEX (L, 1 ) 

MATRIX 

96 

640 

JCOLUM=INDEX(L,2) 

MATRIX 

97 

650 

DO  705  K= 1 , N 

MATRIX 

98 

660 

SWAP=A (K, JROW) 

MATRIX 

99 

670 

A (K,JROW)=A (K,JCOLUM) 

MATRIX 

100 

100 

700 

A  (K, JCOLUM) =SWAP 

MATRIX 

101 

705 

CONTINUE 

MATRIX 

102 

710 

CONTINUE 

MATRIX 

103 

DO  730  K  =  1 , N 

MATRIX 

104 

IF  ( INDEX (K , 3 )  .NE.  1)  GO  TO  715 

MATRIX 

105 

105 

730 

CONTINUE 

MATRIX 

106 

RETURN 

MATRIX 

107 

715 

ID  =  2 

MATRIX 

108 

RETURN 

MATRIX 

109 

END 

MATRIX 

110 
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ANGLE  OP  ATTACK  45*00  DEGREES 


MAXIMUM  RADIUS 
•50000 


OGIVE  LENGTH 
3  •  1  00 


OGIVE  DIMENSIONS 

Z  MATCHING  POINT 
3  ©  I  0  0  0  0 


BODY  LENGTH 

11.00000 


**®**FREE  PARAMETERS***** 

INITIAL  VORTEX  STRENGTH  ©04000 

INITIAL  RADIAL  PERTURBATION  ©00600 

RADIAL  ASYMMETRY  AT  WHICH  FIRST  VORTEX  IS  SHED  ©10000 

PER  CENT  DECREASE  IN  CIRCULATION  OF  SHED  VORTICES  -2*50000 
LOCATION  FACTOR  FOR  NEWLY  INTRODUCED  VORTICES  ©07500 

STROUHAL  NUMBER  ©16660  CORRESPONDING  SHEDDING  DISTANCE  3©0G120 

SEPARATION  ANGLE  CONSTANTS  ©73300  -0*00000  POINT  OF  TRANSITION  1*00000 


INITIAL  STEP  SIZE  ©0002500 


INTEGRATION 
PRINT  FREQUENCY  -0 
ERROR  BOUND 


OPTIONS 

ERROR  BOUND  6 

AFTER  VORTEX  SHEDDING 


PRINT  INCREMENT 
6 


1*0000000 


X  =  *1000000  R  = 

VORTEX  NUMBER 
1 
2 

CL  CML 

©0042325  *0002809 

ERRORS  0©  0© 


0325289  THETAO  =  41©9978X 

LAM  RADIUS 

©0355897  ©0331857 

-©0578202  *0327899 

CN  CMN 

©0000144  ©0000014  0© 

0*  0© 


STEP  SIZE  =  ©25 00000 0E-03 

ANGLE  Y 

46  ©  6493614  ©0241315 

-46© 64936 14  -*0238437 

CD  CY  CYM 

0000000  ©0000029  *0000003 

0* 


Z 

©0227807 

*0225090 


td 

i 


X  =  1*1030758  R  =  ©2956664  THETAO  =  41*99781 

VORTEX  NUMBER  LAM  RADIUS 

1  ©2438641  ©3327822 

2  -©2082865  *3244723 

CL  CML  CN  CMN 

©3496745  ©2415966  ©0430999  ©0372341 

ERRORS  ©607359E-09  ©218640E-09  ©23634QE-Q8  -.6525G9E- 


STEP  SIZE  =  ©21 046273E-0  I 

ANGLE  Y 

39*0463583  ©2096358 

-40*5432500  -©2109141 

CD  CY  CYM 

©1302038  ©0107881  ©0094284 

09  0  © 


Z 

©2584508 

©2465715 


X  = 


ERRORS 


2*1196964  R  =  ©4511470  THETAO  =  41*99781 

VORTEX  NUMBER  LAM  RADIUS 

1  *4714942  ©6214725 

2  -*4356448  ©5754213 

CL  CML  CN  CMN 

8141343  ©9785404  ©4469488  ©7358462 

- ©  665432E-09  -©820530E-09  ©47I271E-10  -©646229E 


STEP  SIZE  «  © 35862 999E-Q 1 

ANGLE  Y 

30  ® 9960132  ©3200449 

-32*8963095  -©3125230 

CD  CY  CYM 

•9504091  ©1361153  ©2284237 

-09  0© 


Z 

*5327282 

©4831553 


X  =  2*8623315  R  =  ©4971352  THETAO  =  41*99781  STEP 

VORTEX  NUMBER  LAM  RADIUS 

1  *6638951  ©8170361 

2  -*6455558  ©7118971 

CL  CML  CN  CMN  CD 

©9885735  1*4009501  1*0517271  2*2520225  0*00000 

ERRORS  - ©  1 05708E-09  *329976E-06  -*291649E-08  -©409414E-08  0* 


SIZE  =  *250  0 0  0  0  0E-03 

ANGLE  Y 

26*8023124  ©3684126 

27.9646625  -©3336277 

CY  CYM 

00  *3855470  .8595407 


Z 

©7292599 

©6287738 


X  =  2*8623315  R  = 

VORTEX  NUMBER 
1 
2 
3 

CL  CML 

©9885735  1*4009501 

ERRORS  0©  0© 


4971352  THETAO 
LAM 

©6472977 

-.6421862 

•0022805 

CN 

1*0344401 

0* 


-  41*99781 

RADIUS 
©8170361 
.7118971 
©5227461 

CMN 

2*2025413 

0* 


STEP  SIZE  =  *  1 0 0 0 0 00 0E-02 

ANGLE  Y 

26*8023124  *3684126 

-27.9646625  -©3336277 

44*9565092  ©3693581 

CD  CY  CYM 

0.0000000  *3598689  ©7860415 

0  *  0  •  0  « 


z 

©7292599 

.6287738 

.3699192 
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X  = 

3.1358214  R  = 

.5000000  THETAO 

=  41.99781 

STEP  SIZE  =  . 17107098E-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

.9293537 

23.6852872 

.3733331 

.8510703 

2 

-.7280791 

.7461542 

-26.2962073 

-.3305552 

.6689390 

3 

.0139121 

•5284685 

46.1529236 

.3811269 

.3660892 

CL  CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000  1.4345633 

1.2155258  \Z 

.7449078 

.9638301  .4855503  1.1628711 

ERRORS  • 166824E-07  .117700E-Q7  -.490723E- 

08  .449925E- 

08  . 344445E-08 

•805208E-08  0. 

X  = 

4.1446835  R  = 

.5000000  THETAO 

=  41 .99781 

STEP  SIZE  =  .35007074E-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

1.4027342 

17.7113434 

.4267421 

1.3362464 

2 

-1.0237597 

.8997440 

-20.3545102 

-.3129559 

.8435625 

3 

.1757541 

.5636118 

39.5839052 

.3591377 

.4343713 

CL  CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000  1.4345633 

1.8985781  5 

.2548438  1 

.3165386  .6803161  1.8155487 

ERRORS  .371962E-08  .239884E-11  .809811E- 

10  -.65<jJ020E- 

09  -.393803E-09  - 

•422120E-09  0. 

X  = 

5.1549628  R  = 

.5000000  THETAO 

=  41.99781 

STEP  SIZE  =  .29757162E-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

1.9788370 

16.2281326 

.5530109 

1.8999934 

2 

-1.2720087 

1.2138898 

-16.4435008 

-.3436154 

1.1642409 

3 

.3921419 

.5724520 

36.5638004 

.3410197 

.4597900 

CL  CML 

CN 

CMN 

CD  CY 

CYM 

1*0000000  1.4345633  2.9112933  9.9895735  1.8015862  -.1686913  -2.2403960 

ERRORS  -.31 1 141E-07  -.123495E-07  .649968E-10  -.286473E-09  -.973207E-12  -.383286E-10  0. 


5.8635320 

R  = 

.5000000  THETAO  = 

41.99781 

STEP  SIZE  =  0. 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

2.4158047 

15.6478473 

.6516006 

2.3262693 

2 

-1.4302811 

1.5175554 

-15.0338253 

-.3936376 

1.4656138 

3 

.5239926 

.5833192 

35.0541973 

.3350300 

.4775104 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

3.7994491  14. 

8912346  0 

.0000000  -1.4771218 

-9.4785667 

ERRORS  -» 165932E-06  -.758859E-07  .640592E-09  .155940E-08  .114145E-08  .486370E-09  0. 


X  a  5.8635320 

R  = 

5000000  THETAO  = 

41.99781 

STEP  SIZE  =  .  10000000E 

-02 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

2.4158047 

15.6478473 

.6516006 

2.3262693 

2 

-1.3945240 

1.5175554 

-15.0338253 

.3936376 

1.4656138 

3 

.5145310 

.5833192 

35.0541973 

.3350300 

.4775104 

4 

-.0028980 

.5263422 

-45.0209574 

.3723162 

.3720440 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

3.7463173  14. 

5796946  0 

.0000000  -1.2954685 

-8.4134367 

ERRORS  0. 

ERRORS  0. 

0. 

0. 

0. 

0.  0. 

0. 

0. 

6.1615320 

R  = 

.5000000  THETAO 

*  41.99781 

STEP  SIZE  =  . 16000  0 OOE 

-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

2.6042565 

15.3967217 

.6914326 

2.5107913 

2 

-1.3945240 

1.7091872 

-13.9904236 

.4132126 

1.6584861 

3 

.5714561 

.5970385 

33.7584075 

.3317697 

.4963707 

4 

-.0076643 

.5335495 

-46.3339277 

.3859569 

.3683915 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

4.0499953  16 

.4051674  1 

.5663684  -1.8966131 

-12.0277866 

ERRORS  .353227E-08  .589238E-09  -.121830E-09  -.504559E-10  .452755E-11  .363392E-11  -.873714E-12  .946037E-12 

ERRORS  0. 


7.1721751 

R  = 

.5000000  THETAO  = 

41.99781 

STEP  SIZE  =  .32094397E 

-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

3.2348758 

14.4416539 

.8067585 

3.1326605 

2 

-1.3945240 

2.4021354 

-11.1526320 

.4646290 

2.3567720 

3 

.8646256 

.7177918 

25.8484516 

.3129517 

.6459770 

4 

-.0825237 

.5600676 

-41.4797306 

.3709636 

.4195971 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

5.0887432  23. 

3413337  1. 

7434051  -3.5828642 

-23.1921245 
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ERRORS  • 132244E-08  -.211717E-08  .292732E-09  -.290729E-10  -.248214E-10  .382313E-10  -.172845E-10  .215497E-10 

ERRORS  0 . 


X  =  8.1991958 

R  = 

.5000000  THETAO  = 

41.99781 

STEP  SIZE  =  .32094397E 

-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

3.8595730 

13.1682765 

.8792562 

3.7580862 

2 

-1.3945240 

3.1316880 

-9.5957769 

.5220401 

3.0878705 

3 

1.1935511 

1.0327161 

18.2391696 

.3232239 

.9808307 

4 

-.2689768 

.5561387 

-38.7058756 

.3477662 

.4339919 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

6.3164263  32. 

7899077  1 

•9861864  -3.7775881 

-24.5257429 

ERRORS  •731932E-09  -.986816E-10  -.325549E-09  -.216827E-09  .257688E-10  .622232E-11  .236255E-11  -.817124E-12 

ERRORS  0. 


X  =  8.8647325 

R  = 

•5000000  THETAO  =  41.99781 

STEP  SIZE  =  0. 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

•6435889 

4.2580371 

12.1574910 

.8967408 

4.1625396 

2 

-1.3945240 

3.6060791 

-9.1234539 

-.5717880 

3.5604585 

3 

1.3749567 

1.3339893 

15.6128605 

.3590245 

1.2847680 

4 

-.3515362 

.5549913 

-38.6311965 

-.3464839 

.4335485 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

7.1938063 

40.2793319  0 

.0000000  -2.9521614 

-17.4474940 

ERRORS  .888596E 

-06  -o  352818E-06  -.744009E-Q9  -.187193E- 

08  .596501E-11  -.720727E-10  -.225967E-10  -.219927E 

ERRORS  0. 

X  =  8.8647325 

R  = 

.5000000  THETAO  =  41.99781 

STEP  SIZE  =  • 10000000E-02 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

4.2580371 

12.1574910 

.8967408 

4.1625396 

2 

-1.3945240 

3.6060791 

-9.1234539 

-.5717880 

3.5604585 

3 

1.3405828 

1.3339893 

15.6128605 

.3590245 

1.2847680 

4 

-.3449832 

.5549913 

-38.6311965 

-.3464839 

.4335485 

5 

.0030946 

.5254779 

44.9261199 

.3710896 

.3720478 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000 

1.4345633 

7.1501010 

39.8918962  0 

.0000000  -3.1014182 

-18.7706156 

ERRORS  0. 

0. 

0. 

0. 

0.  0. 

0. 

0. 

ERRORS  0. 

0. 

0. 

X  =  9.2130058 

R  = 

.5000000  THETAO  =  41.99781 

STEP  SIZE  =  .14707534E-01 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

4.4654066 

11.5864241 

.8968582 

4.3744144 

2 

-1.3945240 

3.8539932 

-8.9834568 

-.6017983 

3.8067181 

3 

1.3405826 

1.5724424 

13.6477708 

.3710216 

1.5280439 

4 

-.3858415 

.5614540 

-37.9124832 

-.3449694 

.4429592 

5 

.0106460 

.5376864 

46.1574387 

.3878041 

.3724441 

CL 

CML 

CN 

CMN 

CD  CY 

CYM 

1.0000000  1.4345633  7.4810645  42.8818751  1.3880138  -2.5306527  -13.6101726 

ERRORS  - ® 165317E-08  .131929E-08  .843059E-11  -.463477E-11  -.817403E-11  .558131E-11  .408816E-12  -.119777E-12 

ERRORS  .639488E-13  -.178651E-12  0. 


X  =  10.2611996  R  =  .5000000  THETAO  =  41.99781 


VORTEX  NUMBER 

LAM 

RADIUS 

1 

.6435889 

5.0895694 

2 

-1.3945240 

4.5943797 

3 

1.3405828 

2.3492428 

4 

-.6270742 

.6291921 

5 

.0965087 

.6233788 

CL 

CML 

CN 

CMN 

1.0000000 

1.4345633 

8.3707261 

51.5497577 

ERRORS  -.141296E-07  -.659408E-08  -.998941E-08  .161486E 

ERRORS  .377907E-10  -.231874E-09  0. 


STEP  SIZE  =  .63391534E 

-01 

ANGLE 

Y 

Z 

9.7266230 

.8598694 

5.0164072 

-8.8936108 

.7102921 

4.5391420 

9.1053566 

.3717686 

2.3196400 

-31.5159190 

.3289010 

.5363831 

40.7982071 

.4073138 

.4719074 

CD  CY 

CYM 

1.4464884  -.7860946 

3.3568254 

-07  -.408549E-09  .977546E-09  .201482E 

>09  - . 208869E-09 

11.0000000  R  = 

.5000000  THETAO  = 

41.99851 

STEP  SIZE  =  0. 

VORTEX  NUMBER 

LAM 

RADIUS 

ANGLE 

Y 

Z 

1 

.6435889 

5.5301893 

8.3632984 

.8043622 

5.4713797 

2 

-1.3945240 

5.1088224 

-8.9891053 

-.7982364 

5.0460762 

3 

1.3405828 

2.9151405 

7.0659801 

.3585980 

2.8930005 
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4  -•8784763  .7550590  -24.7394574  -.3159866  .6857598 

5  .2409266  .6908357  36.0218935  .4062766  .5587426 

CL  CML  CM  CMN  CD  CY  CYM 

1.0000000  1.4345633  9.1400144  59.7364535  0.0000000  .0524773  12.2374875 

ERRORS  . 6444 1 9E-07  .419243E-06  .855520E-07  .533789E-06  .806359E-11  .303639E-11  .688279E-12  .253901E-11 


ERRORS  . 1 60275E-1 1  .12Q698E-11  .892321E-06 


*****PXNAL  STATISTICS***** 

NORMAL  FORCE  COEFFICIENT=CL+CN=  10.14001 
PITCHING  MOMENT  COEFFICIENT=CML*CMN=  61.17102 
CENTER  OF  NORMAL  FORCE=  6.03264 
YAW  FORCE  COEFFICIENT21  .05248 
YAW  MOMENT  COEFFICIENT  =  12.23749 

CENTER  OF  YAW  FORCE=  233.19566 


w 

1 

u> 
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APPENDIX  C 

EXPERIMENTAL  INVESTIGATION 

The  experiments  referred  to  in  this  report  were  carried  out  in 
the  NOL  Supersonic  Tunnel  No.  1.  This  facility  is  an  open  jet, 
blowdown  tunnel,  16  inches  square,  and  is  capable  of  subsonic  as  well 
as  supersonic  speeds. 

The  test  model  had  an  overall  fineness  of  12,  a  nose  fineness  of 
4  and  a  diameter  of  1.1  inches.  It  was  tested  at  the  Mach  numbers  of 
0.5,  0.7,  0.9  and  1.1,  which  corresponded  to  a  Reynolds  number 
variation  of  2(10^)  to  5(10^)  based  on  the  body  diameter.  The  model 
was  swept  from  0  to  70  degrees  angle  of  attack.  Four  runs  were  made 
at  each  Mach  number  with  the  nose  roll  orientation  being  varied  90 
degrees  from  run  to  run.  Select  runs  were  made  using  a  boundary- 
layer  trip.  This  consisted  of  a  grit  strip  1/8-inch  wide  along  the 
leeward  meridian  of  the  body.  Since  the  wind  tunnel  is  an  open  jet 
type,  it  was  possible  to  use  a  subsonic  nozzle  throughout  the  tests 
without  encountering  problems  with  reflected  shock  waves. 

Force  measurements  were  made  using  a  four-component  strain-gage 
balance.  This  type  of  information  was  augmented  by  schlieren 
photographs  which  were  taken  at  increments  of  five  degrees  in  angle  of 
attack.  Oil-flow  studies  were  made  at  each  Mach  number  at  angles  of 
attack  of  20,  35,  50  and  70  degrees. 

Results  from  the  force  tests  indicate  that  normal  force  increases 
slightly  with  increasing  Mach  number,  while  the  yaw  force  decreased 
strongly  under  these  conditions.  The  measured  forces  and  moments  agree 
well  with  the  data  taken  by  Pick  (Ref.  (20))  using  the  same  model. 

Pick's  results,  though,  only  go  up  to  an  angle  of  attack  of  45  degrees. 
The  Strouhal  number  measured  from  the  schlieren  photographs  appears 
to  be  lower  than  that  measured  by  Pick,  although  there  is  considerable 
scatter  in  the  data  obtained  in  the  present  study.  The  boundary- 
layer  trip  significantly  lowered  the  normal  force  and  greatly  reduced 
the  yaw  force  at  the  lower  Mach  numbers. 

In  order  for  force  measurements  to  be  meaningful  at  low  crossflow 

Mach  numbers,  it  is  necessary  to  know  whether  the  boundary  layer  is 

laminar  or  turbulent.  It  is  well  known  that  for  M  $  0.5,  the 

c 

transition  of  the  boundary  is  accompanied  by  a  large  decrease  in  normal 
force.  This  results  from  an  increase  in  the  crossflow  separation 
angle.  At  larger  crossflow  Mach  numbers  the  transition  of  the 
boundary  layer  no  longer  has  a  large  effect  on  normal  force  (Ref.  (3)). 

In  the  tests  conducted,  the  crossflow  Reynolds  number  (Re  sina) 
were  of  the  order  of  10^.  Since  transition  is  expected  near  this  value, 
it  is  difficult  to  predict  the  state  of  the  boundary  layer.  However, 
measurements  of  the  crossflow  separation  angle  taken  from  the  oil-flow 
studies  indicate  that  the  boundary  layer  was  laminar.  This  conclusion 
is  supported  by  the  reduction  of  force  observed  when  the  boundary-layer 
trip  was  used. 
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